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1. Introduction. The standard forms for the general solution of the Confluent Hyper- 
geometric Equation prove too unwieldy for application to many physical problems, 
particularly in the field of Quantum Mechanics. The two standard power series solu- 
tions, .V/,,.,.(y), reduce to a single regular polynomial solution whenever 2m is an 
integer (the standard case for quantum mechanics), and in this case the two integral 
solutions, HW.,.,, (4y), must be computed with an asymptotic expansion which is 


cumbersome for most physically interesting values of y. 

The utility of all of these solutions is further limited because their form necessitates 
undertaking a complete recomputation for every physically significant value of n. This 
makes the labor of computation almost prohibitive in the physically important case 
where both n and y must be treated as continuous variables. 

The possibility of achieving a more manageable form of the solutions was first 
indicated by the work of Wannier® and Jastrow.* Wannier showed that in theory the 
function .V/,,,,(y) could be developed as a series in descending powers of n with coeffi- 
cients given in terms of Bessel functions. Jastrow actually exhibited the first two terms 
of an asymptotically similar series for the solution W,,,,,(y). This paper completes the 
above treatments by producing analytically a general solution of the differential equation 
as a power series in 1/n” with coefficients readily calculable in terms of known functions. 
This treatment differs from those noted above not only in the generality of its results, 
but also in the ease with which successive terms of the series may be explicitly generated. 
The method employed here makes it possible to exhibit the two particular solutions of 
the equation which go to zero as y > 0 and as y > + and to relate these analytically 
to the earlier solutions VW,,,,,(y) and M,,,.(y). Finally this paper will exhibit analytic and 
numerical values for the coefficients of several of the series of greatest physical interest. 

2. The general series solution. For physical applications Whittaker’s standard form 
of the Confluent Hypergeometric Equation, 
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dU 1 n 1/4 — m’ ]., 
—- —_= : Si = ( 
dy + | 4 + - + , it 0, (1) 


is conveniently transformed by the substitutions 


r= 5 ny; m=I1+ ; (2) 
to the form 
- r(t,n) 2 ob a 
c—+ | -4 +2- ED |, 1) = 0, (3) 
dr nN ¥ 2 


is just the hydrogenic radial wave equation in Rydberg units. 


which, with « = —1/n’, 
The form (3) will be taken as standard in this paper. 


The further substitutions 


12 , 1 — 
Z= (8r)'”: ro lea 4 
Br"; iL 9 (4) 
reduce (3) to the form 

(Len ee) i ee e 
Vivo" —n ‘(22) V'™ = 0, (5 

where V, is the Bessel operator of index 2/ + 1, i.e. 
d 5 i . 
V; B2 wettgte — 2+). (6) 


It will be assumed here and proven in the Appendix that the general solution of (5) 
may be written in the form of a power series in 1/n’, i.e. 
7(l,n) —2k yr(l)/ - 
J (2) = dn Vi (2), (7) 
k=0 


where the functions V;’(z) are analytic functions of z and the series (7) converges 
absolutely and uniformly for all real J and for all real n in the region | n | > mo , No being 
an arbitrary positive number. Corresponding solutions of (3) may then be written in 


the form 


Um?) = Don *U;i{), (8) 
in which the U{” (z) are given in terms of the V{’’(z) by the equation 
“(t beta) 
Ui?@ =52Vi ©). (9) 


Since (7) converges uniformly and absolutely it may be inserted in the differential 
equation (5), and differentiated term by term. The terms may then be rearranged and 
the coefficients of the various powers of 1/n* equated to zero, which is necessary and 
sufficient to make (7) a solution of (5). This procedure yields an infinite set of simul- 
taneous differential equations fur the coefficients V;” (z): 

V.iV5’ =0 (10a) 


ViVi? = ( z) Vi, (10b) 
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Now (10a) is just Bessel’s equation of index 21 + 1 whose general solution is the 
cylindrical functions @,;,,(z), an arbitrary linear combination of the Bessel function 
Jo14,(2), and the Weber function, Y,.;,,(z). These cylindrical functions obey the usual 


recursion formulas for Bessel functions: 


2n 
e,,_;(s) + C,.+:(8) _ 2 €,,(z) 
(11) 
z2@/(z) + n@,(z) = 2@,_, (2). 
By utilizing these relations and the definition (8) of the Bessel operator V, , it can be 
shown straightforwardly that 
f2l+2 + q (2 2) d [= 7 (3 " } uf (3 :) 
a! 4(2 re q) Cors+a+¢ 4(3 4. q) 9” Corss+a ia 9° Corsise . (12) 
This equation permits solutions of (10b) to be generated directly. For set g = 0 in 
(12). Then this equation becomes identical with (10b) in the case k = 1, provided that 
V; ’(z) is defined by the bracket on the left side of (12), i.e. by 


wa thet ¥ Fie y 
J 1 = — (4) Cor+3 — 12 (4.) Corse - (13) 





Thus a solution has been found for the first of equations (10b). 

For k = 2, the right hand side of (10b) is just (3z)*V, which, by (13), must contain 
two terms of the form (4z)***@2;.,,, With g = 2 and gq = 3. These are both of the form 
found on the right of equation (12), so that V;"(z) can be generated by two applications 


of the procedure outlined above. Higher terms are generated successively by the same 





process. For example, the first four coefficients of the series for U“’” (z) are 
») 1 
ri ) = (2 z)e, 
1/1 \° 1/1 \* 
* (32) @ - G4) % 
(14) 

oe JAY 11%, L(t.) 

= 16 (32) Cs 30 (3 :) Co + 288 \2~ C; 

« 2 RY : 71 ( :) . 1 (3 :) i ail (! :)" 

= + (32, ©: — § 720 \27) © + 5.760 \27/ © — 10,368 \27/ Se- 


It has been assumed, and will be proven in the Appendix, that the series generated 
above and illustrated in (14) is itself analytic and uniformly convergent in z so that it 
may be differentiated term by term with respect to z to yield a series for the derivatives 
of the solutions of equation (3). Such a series is discussed further in section 5. 

3. Some important particular solutions. The series generated above yields a general 
solution of equation (3) because the linear combination aJ,,(z) + 8Y,,(z) to be inserted 
for the cylindrical function @,, remains entirely arbitrary, so that the coefficients U;’ (z) 
are in fact ambiguously defined. For physical application of the series it is necessary to 
examine the effect of removing this ambiguity by particular choices of the constants 
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oe . P F — — ° ‘* r(l) 
a and 8. This examination is facilitated by defining two sets of coefficients, °U;’ (z) 


r(l ° ° . ° = r(L) ° ° 
and 'U;’(z), which are gained from the ambiguous coefficients U;'’(z) by substituting 
J,,(z) and Y,,(z) respectively for the cylindrical functions @,,(z). Two particular inde- 
pendent solutions of (3), °U“''" (z) and *U (z), may then be defined as the result of 


° . “or r(l lyri) ° 
applying the summation (8) to the coefficients “U;,"’(z) and ‘l : (2), respectively. 
The most general solution of (3) may now be written in the form 


4 7 l ] 4 r — 2h . Orr(l)s. r(l) ~) 
f(l, n) vl (z) + g(l, n) ‘Ll (z) = > n-[f(l, n)°UL? (2) + g(l, n) ‘Uz,’ (2)], (15) 
where f(/, ») and g(l, n) are entirely arbitrary. Important particular solutions of (3) are 
gained by specifying these arbitrary functions in (15). 

Since the Bessel functions all have zeros and the Weber functions all have poles at 


the origin, it can be shown that °U‘'’"’(z) is the only particular solution of (3) with a 
zero at the origin. It must therefore be identical, except in amplitude, with the particular 
solution /,,,.(y) of (1), and a comparison of the leading terms (in z) of the expansions 


of the two series 5 ields 
My,1+12(2"/4n) = no! T(21 + 2)°U"@). (16) 


A second solution of physical interest, the only particular solution which goes to 
zero as z goes to infinity, may be discovered by a comparison of the series developed 


above with a series solution given by Wannier. In the paper previously noted, Wannier 
defines two solutions of (5), J,.,(z) and N2,.,(z), by the formulas 
, a Sy 
J51512) = = M, (2° /4n) 
a Cans 
(17) 
i | [Ta+1+1) ee 
ee) == J51+;(2) cos (21 + 1l)x — J 5,_,(2) |. 
sin (2/ + 1) Lr(n — l)n 


These solutions are shown to be independent and well defined for all values of J and n. 
Wannier further proves that the particular solution which goes to zero at infinity may 


be written 


W. r41/2(2"/4n 2*/4n)? (Pin + 1 + 1)? J2,.,(2) cos (n — 1 — 1) 
2 21+1 
(18) 
+ T(n — In’ Nb14:(2) sin (n — 1 — 1)r]. 
Since the series expansion of M,,.;,;/2 is well known, equations (17) and (18) com- 


pletely determine the expansion of W,,,;.1/2 for any value of 1. It follows that if, for a 
given fixed value of I, the first m + 1 coefficients °U;" and ‘U;” have been developed 
by the generating procedure, the functions f(/, ») and g(l, n) may be determined (to 
terms in ) by explicit comparison of the series expansions (in z) of (15) and (18). 
If 27 + 1 is not an integer, it is convenient to compare the coefficients of the terms in 
in (4z)"'*? in the two series; if 2/ + 1 is an integer the coefficient of the terms 
*!*? log (4z) are most conveniently compared. In the latter case, 2/ + 1 


4z)~?! and 
12)° ** and (32) 
an integer, these two terms of (18) are given by 


5Wannier’s paper has —I'(n — 1)n!*1/2 «++ , a discrepancy which I assume to be due to a misprint in 


the original. 
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; - Tin +1+ 1) (1/22)"" J 
Wa.rsrjole /4n) = +++ + —— (i+ D! cos (n — 1 — 1)r 
ao in (7 l l E lo (3 ) +2 +» + V(n — l) log (n) (19) 
Ss — = ) “a y Z “a — 3 — —_ y e 
7 yn) 6 2 . sam p 





a 21 — r)! 
+ I(n — J(21+ 1)! } i (—1)’ titres es r) = _.} fee, 
r=0 2° (2+1—-rnfa+l+i1—-nr) 
in which ¥(x) = d/dzx log T(x) and y = Euler’s constant. 

This procedure has been carried out for the two most important cases, 1 = 0 and 

1 = 1. In both cases the manipulation yields (to terms in n~"*) 

Wo. ts1o(2?/4n) = n' "Tin + 1 + D[cos (n — 1 — 1)r°U"' 2) 
(20) 
+ sin (n — 1 — 1)r'U"'”(2)]. 


There are additional theoretical reasons for supposing the equation (20) is in fact 
valid for all values of 1, integral and non-integral, but a general proof of this result has 
not yet been given. Until such a proof is produced the method outlined above may be 
used to produce an equivalent result for any value of l for which the coefficients a 
and ‘U;”’ have been developed. 

4. An alternate generating procedure. Since the Weber functions have not been tabu- 
lated for large indices, it is convenient to develop the formulas for Uj’ (z) so that they 
involve only @,(z) and @,(z). This may be accomplished by repeated application of the 
first of the recursion formulas (11) to the functions V{"?(z) generated by the method of 
section 2, but this reduction is arduous and may conveniently be replaced by the pro- 
cedure sketched below. 

V5 (z) is just @2,,,(z), and this may, by application of the recursion 


The function 
formulas, be rewritten in the form 


M 1 2i N (1 2i+1 
Ve (2) = b> a,(4 :) eC, + p bia :) C:, (21) 


‘=m =n 


where the constants a; and b; are known rational numbers and the constants m, n, M, 
and WN are known integers. 
The function V{'’(z) must be expressible in the form 
M+2 1 2i N+2 1 2i+1 
ae ee ae on 


where the constants a; and 8; are unknown rational numbers which can be determined 
by applying the differential equation (10b) to (21) and (22). This application of the 
differential equation is facilitated by the use of the equations 


of ()e] 


vi (3 :)'e, | i me 2) 'e, + iq — 1° — (28 + vid z) @, ’ 


“= 


q@ 


atl = = 1 
~4q(3 :) e, + [¢g — (2l+ y(t :) Cy 


I 


(23) 
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which are immediate consequences of (6) and (11). This procedure reduces V;"?(z) to 
the form (21) after which the method may be reapplied to the generation of V."’(z), ete. 

It must, however, be noted that this procedure, in contrast to that described in 
section 2, does not always uniquely determine the functions V{? (2), for the differential 
equation (10b) is normally satisfied by functions of the form (21) and (22) for all values 
of certain of the coefficients a; and 8; . This will be understood when it is observed 
that if, with V;2\(z) given, a function V;{’(z) is found to satisfy (10b), then the new 
function V;(z) + aV5" (z) will also satisfy (10b) for any a, a fact which follows directly 
from (10a). This ambiguity does not affect the legitimacy of the generating procedure, 
for although the quantity aV;’(z) may be added to any V;’(z), the quantity thus added 
will, by (10b), affect the formulas for V{.,(z), Vi2.(z), ete. It is in fact readily seen 
from (10b) that the net effect of adding aV;'’(z) to the coefficient V;"’(z) is just to 
increase the amplitude of the sum of the series, i.e. of V“'"" (z), by the factor (1 + a/n”). 

This new generating procedure was used in the preparation of the tables which 
follow. For simplicity of computation all those coefficients, a; and 8; , which were not 
explicitly determined by the procedure were set equal to zero, thus reducing the com- 
plexity of the formulas for the coefficients V;;"’(z). It is readily seen that, for 1 = 0, the 
tabulated coefficients are just those which would have been gained using the original 
generating procedure. For / = 1 the tabulated coefficients differ from those provided by 
the original procedure, but in this case the sums of the series gained with the tabulated 
coefficients may be made equal to the sums of the series gained with the standard co- 
efficients described in section 2 by multiplying the former with the amplitude factor 
n’/(n®> — 1). This modification of the amplitudes of the solutions is of no significance 
except when it is necessary to use (16) and (20) for explicit computation of Whittaker’s 
functions, VM, ,, and W,.» . 

5. The tables. Table I below lists the formulas for U, 
Table II gives a similar list for the functions U;” (z). Tables III and IV list corresponding 


formulas for the functions D,"’(z) and D,’(z) which are defined by 


"(z) with k = 0 through 7. 


( 


l — ] ; d r(l)s ‘ 
D;, '() = (2 :) i U, (2), (24) 
from which the functions 
, - Seay ii 
D2) = Don Di = ; oe (25) 
k=0 az 


may be computed. 

Tables V through VIII list values of the functions °U;"’(z), "Uz? (z), °Dz°(z), and 
"D;” (z) for 1 = O and 1, k = 0 through 7, and z = 3.5(0.5)7.5. As before a superscript 
zero preceding the function indicates substitution of a Bessel function for the corre- 
sponding cylindrical function, and a superscript one connotes use of the Weber function. 
Persons interested in utilizing these numerical results may also find useful the Tables of 
Coulomb Wave Functions recently prepared by the National Applied Mathematics 
Laboratory of the National Bureau of Standards. These supply a single irregular solution 
of the wave equation for 1 = 0 and for values of z smaller than those listed in our tables. 

6. Acknowledgments. The problem treated in this paper arose during the preparation 
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Since completing this paper, my attention has been called to the fact that the gene- 
rating procedure developed in section 2 above has been previously discussed by Yost, 
Wheeler, and Breit® in treating the formally identical problem of the wave equation 
with a repulsive Coulomb field (proton-proton scattering). The additional analytic and 
numerical material displayed here seems more than sufficient to justify the presentation 
of this independent investigation. 

APPENDIX 
. . . 2 r . . nA ¢ 
The existence of an expansion in 1/n°. We may seek a solution of equation (5) in the 
form of a power series in 2, 
r(l,n) 
ve"'(@) = 2" Sas (26) 
Gut 
Throughout its circle of convergence this series may be inserted in the differential equa- 
tion, and manipulations identical with those used in producing solutions of Bessel’s 
equation show that the two series 


i 1 #(21+1) @ : 1 2k “ 
(z) = (3 :) » (-1) a,(3 :) (27) 


k=0 


are solutions of (5) for an arbitrary value of a) , providing that the coefficients a, are 


generated by the formulas 


a, = a[1 + (21 + 1)]”’ (28a) 


«, = S-t- Seal) (28b) 


kfk + (21+ 1)) 
In equations (28a) and (28b) the positive or the negative signs are to be taken together, 


* 2] + 1 is not an integer, these equations define two independent solutions of 


cothas 





so that, 1 
(5). That cheice of sign which makes +(2/ + 1) > 0 produces a solution which is regular 
0; the other choice of sign yields a solution which is irregular there. These two 


at 
“regular” and the “‘irregular’’ series, 


solutions will hereafter be distinguished as the 
respect ively. 

When 2] + 1 is an integer, equations (27) and (28) define only the regular solution 
of the equation, but in this event a second irregular solution may be defined by any of 
the usual’ devices developed for Bessel functions. In this paper we make the expedient 
assumption that 2] + fad is not an integer and then produce two independent solutions 
of (5), °V“"''(z) and ' (z), which remain finite and independent as 2/ + 1 is varied 
continual th aia any integer or zero. Because of their continuity as functions of / 
these solutions may be assumed valid for 2/ + 1 an integer, and this assumption may 
be rigorously justified by standard methods. 

We now investigate the convergence of the two series subject to the simplifying re- 
striction that m and 1 be real variables, z remaining complex. It then appears from (28b) 
that all the coefficients a, of the regular series have the same sign as a» , so it follows 
directly that 





6F. L. Yost, J. A. Wheeler, and G. Breit, Phys. Rev. 49, p. 174 (1936). 
7G. N. Watson, Theory of Bessel functions, (Cambridge University Press, 1945), second edition, 
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8 
e, i hfe + (21 + 1)] 
—— =: 2 Cee Brae Re: <r ase si EE 
a  (k+)ikK+1+ 21+] * 4? + (a2/u-) 
(29) 
l l 
<7 LIL aa. cat 2 
(k+ Dik +1+ (21+) ° 4n 
The inequality in (29) provides a very weak condition on a,,,/a, , because it is gained 
by dropping the entire positive term a,-2/4n’a,_, from the positive denominator 
1 + a,_./4n*a,_; in the equality. Stronger conditions may be gained by successive appli- 


sation of (29) to itself, for (29) may be used to set a lower bound on a,_,/a,-, , and this 
bound may be used to write a second and stronger form of (29). The process may be 
repeated indefinitely, and the p-th form of (29) gained in this manner is 


a _..._ (g9+%) Poe om i 30) 
a, 2k +1 — 2p)[kK +1 — 2p 2 (21+ 1] (p + 1)4n° 
a formula which reduces to (29) when p = O and which is valid for all k > 
2p — 1 = (21 + 1) 


It follows that for all || > ny and for all |z| < R (where n, and R are arbitrary 
positive numbers) the regular series in (27) converges absolutely and uniformly, for 
if p in (30) is chosen to be the largest integer less than R°/16n, , an integer ky may always 
be found such that a,.,; | $2 \"/a, < 1 for all k > k, . It follows that the regular series 
defined by (27) converges to an analytic function of z in any bounded domain in the z 
plane and that this function is uniformly continuous in n for all real n such that | n | > 
fy > VU. 

It may next be noted that equation (28) generates coefficients a, which are given by 


finite polynomials of the form 
Vi 
k —2m ‘ 
a, = >, bin 31) 


where M/, = 4k or 4(k — 1), whichever is an integer. Further, if a, is a coefficient of 
the regular series, the polynomial coefficients b,, must be positive quantities. 

3y considering (27) and (31), we may now define the ‘‘complete series” for V (z) 
as the series in which each b,, n-*” (4z)**~"*'*” is considered a separate term and in 
which for each value of k, the summation is carried out over values of m from 0 to M, 
before k is increased by unity. It then follows from the absolute convergence of (27) 
and from the uniformly positive values of the b;, that the “complete series” for the 
regular solution also converges absolutely and uniformly, so that the terms of the 
“complete” series may be rearranged to provide an analytic solution of (5) in the form 


VO") = )> Vi?(@)n-™* (7) 


The functions V;"’(z) are defined by uniformly and absolutely convergent power series 
in z, so that they are analytic functions of z, and the entire series (7) converges uniformly 
whenever || > %%. 

The existence of the series (7) has so far been proven only in the case of the regular 
solution of (5), but the proof is readily extended to the irregular solution and hence to 
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the general solution, which is an arbitrary linear combination of the two. For consider 
a series (27) in which the coefficients a, are replaced by coefficients aj determined by 


at = af | 1 + (21 +1) [7 (32a) 





they + (ab/4n) 
uF kt Ql+)| (32b) 


This series contains only positive terms, so that the entire proof applied to the regular 
series holds for it. Further, if the choice of sign which produced the irregular series in 
(27) is utilized in (32), then every aj and every b/* is greater than or equal to the corre- 
sponding a, or b%, of the irregular series. It follows from the comparison test that the 
“complete series” for the irregular solution must converge absolutely and uniformly, so 
that the rearrangement of terms is again permissible. This is sufficient to justify entirely 
the assumption made in section 2, above. 


Taste I. Formulas for the coefficients U;°’(z) with k = 0 through 7. 


. 2 2 4 1 1 2 2 4 3 
(2) all Las @) + a (2) lee C1 
it. 2 Lig.  m iee 1 (2) 
“= —| 555 ~ 252 (3) + 60,480 (2) - 10, sa ( (2) (2) Co + E 168 \2 
» 18 (st = 13 me |G ‘Ve, 
' 60,480 \2 25, 920 \2 2 
r>- [+ a) ~ st) © wT) 
‘ “™ 1240 240 \2 12,096 \2 saa sa 2 rh ee 
-[ 35-8) +35 - ae) + (3) 
240 +160 \2 12 \2 at + Tees 2 
7 497, ai ( (2)" (ye. 
~— 1 (z\ z= i§ (2) 1,977 = (2) 
Us 4-3 ~ 732 0) +7 21, (3) - 3,448 \2/ + 31,933,440 \2 
- waitan(®) + waa G10) 
43, 545,600 2 T 39 859,840 2 a ™ 
Lia a+ ala 
132 +~=«88 \2 63, 360 10,560 \2 


8 
479 (zt. 701. (2)" + 13 (2)"|(Ve 
2,128,896 \2 43,545,600 \2 29,859,840 \2 y7™ 


Cin 
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yo =|; 691 691 ay «sm (2) _ _ 3,313 (2) 43,037 (2) 
é 32.760 32,760 \2 27,027 \2 1,853,280 \2/ * 197,683,200 








_ 


160,361 ay 509 (ey 1 (2)""| 4} 
= ane « 7 = Pergamon Co 
10. 378 368, 000 Z 870,912 000 2 119,439,360 \2 \2, 

; 691  —_—<691_— (s) + 14,929 (2) 6,977 (2) 
32,760 21,840 \2, * 864,864 \2 e 182 624 
ae 47, 951 (2) -. so) + 461, 819 (2) 
65, 894. 100 \2 5.930 196,000 \2 134,120,448,000 \2 


7 aes ca 600 (2) + arr 908, 480 (2): ‘| Zy'e, 





10) 273 273 4) 26,609 (2) 6,953 (2) 3,569 (2) 
~~ ae 3,276 \2 T 7 786,240 943,488 \2 + 3,706, 6.560 \2 ) 
ai (z)" 4 —404,561__ (2)" 7 65,539 (2)" 
148,262,400 \2 99 632,332,800 \2 536,481,.792,000 \2 
11 — I (2)" |(2) ' 
+ 6,270.566,400 z)" ~ 180,592,312,320 9) 2 Co 
. | 273 ue (2) 53,909 (2) 45,011 (2) _ 5,363 (2) 
| 7 - ~ 9 18 34 \2 a, 786,240 \2/ ~ 2,358,720 + 1,729,728 \5 


\14 


. 279, 187 [sy — 286,511 7 (2)" 919,637 _ (2) 
889,574,400 \2 14,233,190,400 \2 1,162,377,216,000 
353 i - 61 (2)"|(2) , 
7 20,901, 888,000 451, 480, 780, 800 \2 2 bie 


TaB.e II. Formulas for the coefficients U;"’(z) with k = 0 through 7. 





(on | 
ll 
jet | 
n|™ 
"ita 
DNln 
ed 
re 
| 
—_ 
dole 
ee 
© 


aa. coe oma) = (2)' (ze 
vse iia = 20,160 \2 + 10.368 \2 2 Co 
ry | 31 i: _ 31 iz) ai = (2)' (Ze 
15,120 12,096 \2/ * 1,440 \2/ j\2/ * 
yo =F . (2) +- 481 oe __ 4B (2) |(z)'e _ fia 
+ ™ | 30,240 40,320 \2 1,814,400 \2 622,080 \2/ |\2/ ~° ~~ 30,240 


_ 4 (2) 799 | (2) 2,111 (2) 1 fz ‘ z) 
24,192 \2 ¥ 1,209,600 2) 21 772,800 2 + 497,664 (5) (3 ws 
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ro “| 31. (2) 4 —25_ (2) — 5:88 5,557 (2) + tt (2) 
_ 15,840 21,120 \2 59,136 \2 95, 95.800, 320 \2 3,225,600 \2 





1 {z\'° |/2z\° 31 31 z\’ 661 Z P 
- saa ()"19'e+ [atts — ata +a 
29,859,840 \2 2/ 15,840 12,672 \2 665,280 2 
_ 3,599 f(z)’, 1471 (zy aotasa (2) Ki)" 
19,160,064 (3) +r 87,091,200 \2 1,990,656 \2 2 C1 
10,331 _ 10,331 (*) +; ret (3 
2,162,160 2.882.880 \2 74,131,200 
24.9% z\* 3,239 - 2z\'"? |/z\° 
+ 2hsl__(2)_ 3.2302)" 4 1 (2)" [Ve 
i.761.177.600 \2 5,225,472,000 \2 107,495,424 \2 2 
=! 10,331 10,331 = (2) + 2,251_ (2) _ _ 147,967 _ (2) 
2,162,160 1,729,728 \2 915,200 \2 296,524,800 \2 


r 
, 6: D0, 857 (2) 229,793 (2) 4 — (2)" 
T 11,623,772,160 \2/ ~ 67,060,224,000 \2 2,351,462,400 \2 


Nin 
“eeaess™ 
es 
AE fo 
i | GO 
lon 
= | 
3/9 
S| 
a, 
bole 
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~ 
| 
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bo 











2,1 my 908,480 (3) |e 
ra. _| 3.421 _ _3,421 (2) + 40,907 _ (2) _ _ 57,131 __ (2) 
7 ~ ~L 196,560 262,080 \2 10,378,368 \2 88,957,440 \2 
1 2,748 (2) - 1,837,343 (s 4 +r _ 200,159 (2)" 
13,929,600 \2 498, 161,664,000 109,445,376,000 \2 


Cin 


(ae +| 28h 

522 2 196,560 

_ 3421 (2') , _ 9,749 (3 235,111 (z)" + 819,631 _ (2) 
157,248 2) 1,081,080 124,540,416 \2 3,558,297,600 \2 

_ __ 407,413 (2)" 4 800,593 _ (2)" _. _ 7a = (2) 
23,721,984,000 \2 1,046, 139,494,400 \2 41,803,776,000 \2 

- asatacpn (1 

" 75,246,796,800 \2 i 


Tae III. Formulas for the coefficients D,° (z) with k = 0 through 7. 
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7,911,653 _ (2)" + 4,001 (2)" 
~ 3,487,131,648,000 \2 134,120,448,000 \2 


13 “)" 1 (2)"|(: 3 
+ 39,948,627,200 (3 ~ 180,592,312,320 \2 ;) Cs 
Tas.e IV. Formulas for the coefficients Dj” (z) with k = 0 through 7. 
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4 313,097 __ (z) - 308,233 (z)" _ ___198 (2) 
| 7,264,857,600 \2 201, 180,672,000 \2 18,811,699,200 \2 
e (2) ‘(ae 2\* -| 10,331 10,331. (2) 
| vy} 480 \2 Co 1,081,080 786,240 \2 


15,079 (3) ~3 423,751 _ (2)' + 992,969 = (3) 
7,207,200 \2 296,524,800 \2 5,811,886,080 \2 
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~ \10 w \ 12 
, ” 6. d 718, 464, 000 


,285,159 
43 5,891, 156,000 


1] 2 
7 2,149,908, (9 





Dp” = _| 3,421 3,421 (2) 4 42,671 (2) 1,223,947 8] 
on. fF 112,320 \2 3,991,680 \2 622,702,080 \2 





14,718 (s) __ 558,461 (2)" 7,406,743 (2) 
71 165,952 \2 45,287,424,000 \2/ * 30,922,789, 888,000 


7 (#\" il (2)" (2) 
a ee — =I 
5,971,968 ,000 3) 128 994. 508, 800 \2/ 2 
’ | 3.421 __ 37,631 (2)' 4 39,815 (2) 52,667 (2) 
“ 198,280 786,240 \2/ * 1,729,728 \2/ 9 580,032 \2 
, _1,316,273 (2) 1,370,671 Ay F 25,463,237 (2)" 
' 1,779,148,800 \2 23.721,984,000 \2 " 10,461, 394,944,000 \2 


a 961 (z)" = 17 (2)" 
25,751,126,016 \2/ 37,623,: 398, 400 


1 z\'8 (2)'e 
T 780,592,312,3 = (3) i _ 














TaBLe V. Values of the coefficients °U{" (z) and 'U;,°’(z) with k = 0 through 7 and with 
z = 3.5(0.5)7.5 
0 2 3 4 5 6 7 
3.5 | + .2404 + (3585 1 0593 + 0040 + 0001 
4.0} — .1321 1855 | + .1662 | + .0204 | + .0013 | °U, (2) 
45|— .5199 1653 | + .3677 | + .0803 | + .0088 | + .0006 | - 
5.0} — .8189 | + .1516 | + .6405 | + .2512 | + .0451 | + .0048 | 
5.5 | — .9390 5591 + 8194 + .6326 | + .1844 + .0304 | + .0033 
6.0 | — .8301 1.6394 | + .5094 | +1.2653 | + .6126 | + .1543 | + .0248 | + .0028 
6.5; — .5000 2.8582 | — .9198 | +1.8590 | +1.6560 | + .6394 | + .1500 | + 0241 
7.0| — .0016 | —3.7693 | —4.1461 | +1.2775 | +3.5502 | +2.1826 | + .7462 | + .1703 
7.5 | + .5072 3.7948 | —9.4364 | —3.0898 | +5.4044 | +6.0816 | +3.0829 | + .9871 
3.5 | + .7178 | + .03855 | — .0362 | — .0043 |} — .0003 —_ 
4.0/ + .7959 | + .2879 | — .0374 | — .0188 | — .0020 | 1U,(z) 
4.5 | + .6772| + .7016 | + .0522 | — .0286 | — .0073 | + .0009 — - 
5.0 | + .3697 | +1.1968 | + .3751 | — .0125 | — .0200 ;} — .0018 
5.5 | — .0653 | +1.5768 1.1173 | + .1787 | — .0276 | — .0124 | — .0067 
6.0 | — .5250 | +1.5515 | +2.3930 | + .9167 | + .0838 | — .0297 | — .0128/} + .0091 
6.5 | — .9908 | + .8266 | +4.0212 | +2.8948 | + .7845 | + .0375 | — .0299 | — .0127 
7.0 | —1.0593 | — .7569 | +5.2477 | +6.9946 | +3.4453 | + .7529 | + .0207 | — .0512 
7.5|— ‘ial: —3.0720 | +4.5695 |+13.5834 venue +4.3933 | + sai — .0145 
| | | 


| 
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TABLE VI. Values of the coefficients °U{"(z) and ‘Uj (z) with k = 0 through 7 and with 
z = 3.5(0.5)7.0 

















\Fl 0 1 | 2 | 3 4 5 6 7 
3.5| + .6768 | — .4812| — .1773 | — .0175 | — .0008 

4.0| + .8603 | — .3974| — .3854 | — .0712 | — .0060 | — .0003 U4 (z) 

4.5) + .9556 — .0266 — .6660 — .2273 — .0328 — .0027 ‘eee 
5.0| + .9121 | + .7015 | — .8679 | — .5861 | — .1402 | — .0184 | — .0015 

5.5 | + .7043 | +1.7426 | — .6349 | —1.2141 | — .4851 | — .0990 | — .0126 

6.0 | + .3443 | +2.8845 | + .6039 | —1.9323 | —1.3738 | — .4327 | — .0822| — .0106 
6.5 | — .1149 | +3.7384 | +3.4891 | —1.8745 | —3.1495 | —1.5606 | — .4354 | — .2526 
7.0 | — .5864 | +3.8027 | +8.3606 | +1.0917 | —5.5007 | —4.6443 | —1.9128 | — .5038 








- .0143 


a. 4 











3.5 | — .6271 | — .4019 | + .0517 | + .0024 

4.0| — .3640 | — .8184| — .0342 | + .0343 | + .0067 | — .0019 117, )(z) 

4.5 | — .0203 | —1.2730 | — .3460 | + .0369 | + .0196 | + .0006 : 
5.0 | + .3657 | —1.5819 | —1.0580 | — .0954 | + .0392 | + .0102 | + .0052 

5.5 | + .7270 | —1.4944 | —2.2842 | — .6801 | — .0028 | + .0331 | + .0081 

6.0 | + .9847 | — .7640 | —3.8890 | —2.3339 | — .4337 | + .0339 | + .0256 | + .0134 
6.5 | +1.0686 | + .7416 | —5.2455 | —5.9138 | —2.2843 | — .2992 | + .0439 | + .1615 
7.0 | + .9383 | +2.9197 | —5.0617 oe —8.0027 | —2.3942 | — .2510 | + .0885 











TABLE VII. Values of the coefficients °D{° (z) and *D,°’(z) with k = 0 through 7 and with 
z = 3.5(0.5)7.5 



































| ) | 
h 0 l 2 3 4 | 5 6 | 7 
i — “apie wat ~ gloat angel -pmee SS ee ee | 
3.5 | —1.1641 | + .5464 | + .2399 | + .0248 | + .0012 ———— 
4.0 — 1.5886 + 3094 + .6046 + .1219 + .0108 + .0005 oP), (z) 
1.5 | —1.6227 | — .6451 | +1.1299 | + .4524 | + .0696 | + .0059 ——_—— 
5.0 | —1.1100 | —2.5143 | +1.3511 | +1.2957 | + .3423 | + .0473 | + .0041 | 
5.5 | — .0518 | —5.0341 | + .1217 | +2.8315 | +1.3213 | + .2896 | + .0387 
6.0 | +1.3558 | —7.2422 | —4.5006 | +4.2666 | +4.0382 | +1.4073 | + .2830 | + .0381 
6.5 | +2.7472 | —7.5067 |—14.4887 | +1.9408 | +9.5432 | +5.5044 | +1.6547 | + .3222 
7.0 | +3.6760 | —3.9742 |—29.9598 |—13.2338 |+15.4978 |+17.2609 | +7.8914 | +2.2065 
7.5 | +3.7454 | +4.5633 |—46.3111 |—56.8456 | +5.9852 | 41.6448 |+30.6623 |+12.3593 
| ee es oe — 
| | | 
3.5 | + .5789 | + .5965 | — .0599 | — .0197 | — .0035 , 
1.0] — .0678 | +1.3490 | + .1059 | — .0547 | — .0115 | + .0030 1D, (z) 
4.5 | — .9857 | +2.1480 | + .7959 | — .0496 | — .0386 | — .0030 —— 
5.0 | —1.9282 | +2.4001 | +2.5091 | + .3255 | — .0799 | — .0280 | — .0042 
5.5 | —2.5673 | +1.2654 | +5.5298 | +2.0747 | + .0874 | — .0904 | — .0119 
6.0 | —2.5938 | —1.9924 | +9.1250 | +7.2941 | +1.6581 | — .0485 | — .0705 | — .0921 
6.5 | —1.8299 | —7.4553 |+10.4872 |+18.8387 | +8.7051 | +1.4192 | — .1036 | — .1935 
7.0 | — .3179 |—14.0041 | +4.1085 |4+37.7230 |+31.3594 |+10.8055 |+1.4318 | — .2809 
ed +1.6497 |—19.0857 |—17.2733 we +87.7617 |+51.1071 nese +1.7063 
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TABLE VILL. Values of the coefficients °D{” (z) and 'D,” (z) with k = o through 7 and with 
z = 3.5(0.5)7.0 











\ | | | l 
Nd ee | Seen oF ce © Ff | 6 | 7 
3.5| + .7277| — .0749 | — .5522 | — .0941 | — .0063 | | ——_——_ 
10 | + 5962 | + 8377] —1.0240 | — .3628 | — .0441 | — .0031] = | DiC) | 
45\+ 1473 | +2.4529 | —1.2668 | —1.0756 | — .2312 | — .0281 | —_— 
5.0| — .6210 | +4.5236 | — .3164 | —2.4620 | — .9419 | — .1797 | — .0172 


5.5 | —1.5915 | +6.2828 | +3.4468 | —4.0832 | —3.0436 | — .9013 | — .1378 
6.0 | —2.5302 | +6.4536 |+11.7834 | —3.3613 | —7.7458 | —3.3690 | — .8705 — .1361 
6.5 | —3.1324 L3 5506 |+25.1749 | +6.6506 |—14.5558 | — 12.1739 | —4.4488 | —1.0093 


wo 


ww 








7.0 | —3.1059 | —3.4974 |+40.5401 |+38.8142 |—14.2714 |—32. 1868 —18.6006 | —6.1044 

3.5 | + .7660 | —1.2584 | — .0823 | + .0565 | + .0074 | | —_—— ) 

4.0 | +1.2277 | —1.8344 | — .6849 | + .0815 | + .0304 | + .0051 | | 1D,(z) 

4.5 | +1.6832 | —1.8877 | —2.1735 | — .1481 | + .0840 | + .0260 Lnetipeicend 

5.0 | +1.9322 | — .7766 | —4.8160 | —1.4141 | + .0598 | + .0853 ze .0010 | 

5.5 | +1.7750 | +2.0529 | —8.0904 | —5.4407 | — .7879 | + .1463 | + .0442 

6.0 | +1.0840 | +6.6706 | —9.8249 |-14.8121 | —5.2657 | — .4507 | + .1660 | + 1452 
| 


6.5 | — .1295 |+12.2448 | —5.5858 | —31.2672 —20.9246 | —5.4457 | .1685 | + .2918 
7.0 | —1.6797 |+16.8159 +10.8370 |—50.8593 |—62.7531 —29.3100 | —5.7721 | + .0160 
| | | 
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THE PRINCIPLE OF MINIMIZED ITERATIONS IN THE SOLUTION OF 
THE MATRIX EIGENVALUE PROBLEM’ 


BY 
W. E. ARNOLDI 
Hamilton Standard Division, United Aircraft Corporation, East Hartford, Conn. 


An interpretation of Dr. Cornelius Lanczos’ iteration method, which he has named 
“minimized iterations”, is discussed in this article, expounding the method as applied 
to the solution of the characteristic matrix equations both in homogeneous and non- 
homogeneous form. This interpretation leads to a variation of the Lanczos procedure 
which may frequently be advantageous by virtue of reducing the volume of numerical 
work in practical applications. Both methods employ essentially the same algorithm, 
requiring the generation of a series of orthogonal functions through which a simple 
matrix equation of reduced order is established. The reduced matrix equation may be 
solved directly in terms of certain polynomial functions obtained in conjunction with 
the generated orthogonal functions, and the convergence of the solution may be observed 
as the order of the reduced matrix is successively increased with the order of the original 
matrix as a limit. The method of minimized iterations is recommended as a rapid means 
for determining a small number of the larger eigenvalues and modal columns of a large 
matrix and as a desirable alternative for various series expansions of the Fredholm 
problem. 

1. The conventional iterative procedures. It is frequently required that real latent 
roots, or eigenvalues, and modal columns be determined for a real numerical matrix, u, 
of order, n, in the characteristic homogeneous equation,* 


QAI — wk =0 (1) 
which is satisfied by any of n values of the scalar, \ = A, with their associated modal 
columns, k = k, . Beginning with an arbitrary column, k,) , and repeatedly premulti- 
plying by u, the iterative procedure, k;,, = uk; , will converge to the modal column 


corresponding to the largest, or dominant, latent root. After obtaining this solution, the 
dominant mode may be removed by any of several methods, so that the next largest 
root of the original matrix becomes the dominant root of an altered matrix, whereupon 
the same procedure may be repeated to obtain the next root, and so on, until all desired 
roots and modal columns have been obtained. Since the accuracy of each root and modal 
column is dependent upon the accuracy with which each previous column has been 
determined, this method obviously requires a wasteful amount of labor if the initial 
roots and modes are not needed, as is often the case in practical applications, and con- 
vergence can be extremely slow if the roots are not widely dispersed. 

The non-homogeneous equation, or Fredholm problem in matrix form, may be 
represented by 


Ql — uk = q, (2) 


1Received May 15, 1950. 

*The general notation used in this discussion to represent matrices, columns, rows, and matrix 
equations follows as closely as possible the conventions established in ‘Elementary Matrices”, Frazer, 
Duncan and Collar, Cambridge, 1938. Brackets and braces are employed only when necessary to avoid 
confusion, e.g., to distinguish between the column, k; , and the matrix of m such columns, [k,]. 
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—_ 
i? 2) 


where \ and g are specified and a solution is desired for k. This solution may be obtained 
by direct treatment of the array as simultaneous equations, a common algorithm being 
the method of pivotal condensation [1], a lengthy and laborious procedure but some- 
times indispensable. Among essentially iterative methods are the Schmidt expansion [2], 
which requires first the solution for latent roots, modal columns, and modal rows, 


whereupon k is found as 


] K,¢ 
k = z. ——— *rd p, (x, = modal row). 

— k= A. Bk, 
Another iterative method is based upon the Liouville-Neumann expansion [2], where 


the reciprocal matrix is written as an infinite geometric series, 
k= (J —u)"¢q =A 10 +N ‘utr w+ ---)dq 
whose solution is obtained by iteration: 
Kies = NX 'g + Nuk; 


Unfortunately, this method converges only if \ is greater in magnitude than the magni- 
tude of the dominant root of u. It can also be applied in special cases when the field of 
roots can be shifted so as to meet the requirement for convergence, but the method is 
nevertheless of limited application and is often subject to slow convergence. 

2. A classical method for reducing matrix order. The solution of the original 
homogeneous equation may be replaced by the solution of a matrix equation of reduced 
order by introducing the approximation that the modal column, k, may be represented 
by the sum of a series of arbitrary columns, k; , multiplied by coefficients, c, , to be 
established in the solution. Substituting 


m 


k= Dek; = [kile (3) 


1=1 

into Eq. (1) yields 

(AI — u)[kiJe = 0, (4) 
where [k;] is a rectangular matrix whose columns are k, , and c is a column of elements 
c; . Now premultiply Eq. (4) by another rectangular matrix whose rows are x; , corre- 
sponding with k; in such fashion that the row x = a o c,k; Will satisfy the original 
matrix in the form, 

K(AL — u) = 0. 


It should be noted that the functions, k; and x; , are also required to satisfy the point 
boundary conditions implicit in the matrix, u. The matrix equation now becomes 


[x,](Al — u)[k,Je = 0, 
or, in more concise form, 
(Al — [k,ki]"[k;uk,])c = 0. (5) 


Equation (5) will henceforth be referred to as the “reduced equation,” and the matrix, 
[x;k;]"'[x;uk,], will be known as the “reduced matrix.” Since m functions were assumed 
to be sufficient for an approximate solution, this equation involves a matrix of m order, 
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presumably less than the original n order, hence more rapidly solved than the original 
by conventional iterative procedures. If one could hope that m could be chosen sufficiently 
small, a solution might be available by direct expansion of the determinant, but it is 
usually necessary, in the general form given above, to employ a relatively large number 
of functions, k; and x; , so that the advantages of this transformation of the original 
matrix are not universally apparent. Furthermore, the labor of calculating the elements 
of the reduced matrix may be great, particularly since the inversion of [x;k,] is required. 

The Galerkin method [3] is a variation of the process outlined above, differing only 
in the choice of x; . By the Galerkin method, each x; is the transposed of the corresponding 
column, k; , which somewhat simplifies the determination of the inverse matrix re- 
quired. However, it is also objectionable on the basis of tedium. 

3. The Lanczos method of minimized iterations. The following exposition differs in 
form from the discussion originally given by Lanezos [2], but it arrives at the same 
results and makes a further extension of the Galerkin type more evident. Lanczos 
reduces the matrix order as described above but eliminates the objections of labor in 
the formation of the reduced matrix, indicates a solution in terms of polynomial equa- 
tions derived from a direct iterative procedure involving the original matrix, avoids the 
necessity for separately forming and solving the reduced matrix, and provides a means 
for efficiently generating the required columns and rows in order that a convergent 
solution will be obtained with a minimum number. 

To deduce the Lanezos method of minimized iterations, it is first observed that the 
formation of the reduced matrix would be greatly simplified if the matrix, [«;k;], were 
of diagonal form, the otherwise arbitrary rows and columns initially chosen being of 
such form as to satisfy x;k; = 0 when 7 # 7. The inverse would then consist of a diagonal 
matrix whose individual elements would be the reciprocals of the corresponding elements 
in [x;k,], and the reduced equation would take the form 


(xz ~ Ea = 0. (6) 


This form has further potential advantages, obtained by observing that the desired 
orthogonality relationship among the chosen rows and columns may be established by 
generating these rows and columns from an initial row and column according to the 


following expressions: 


Kiss = uk; oe ak; — B; _k; — Vi-2ki-2 vecd b;<3k;~3 Som ee ae (7) 
ies = ci— A;K; sine B;-1«; ie Ci-2Ki-2 ‘aac D;~3ki-3 mae ya (8) 
The necessary scalar constants, a; , Bi-1 , Yi-2, *** » Ai, Bi-a , Ci-2 , +++ , may be de- 
fined by forming certain scalars as row-column products and noting that many elements 
vanish, due to the orthogonality relation postulated, «,k; = 0 for i # j, as follows: 
Kikea. = xuk; — a,x,k; = 0, 


Kiki 42 _ KU; 44 ou Bi xk; = 0, 
Kikias aad KUK; 49 — vikik; - 0, 


etc. 








Note that a, 


uk; 


Keele, 


KxU 


Kuk; +1 


Therefore, 


B; 
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also 
Kiwik; = wuk; — A;«,k; = 0, 
Kiaok; = Kuk; — B;xk; = 0, 
Ki+3ks = Ki42uk; — Ci Kxsky = 0, 
etc 
whence 
Kuk, xuk; 7) 
a=-—, A, = — y=, etc., 
Kk, Kk; K,h; 
uk, Kea Uk, K;40uk 
A, = B, = =i C, = “Ges ete. 
k,k; ik, Rik, ° 


A, . Other identities may be found as follows: 


Riss + a;K; + Bi-1ks-1 + ar ae 


Kia kiss ° 


Kier tb Aces + Berk t+ -::, 


Keatkiay . 


= B,. 


The later constants will vanish, as is next shown. 


cu = 
Therefore, 
Kukseg = 0 
whence, 
y: = 0. 


Similarly, C; 
vanish. 


The functional relationships by 


0, and all further scalars, 6; , €, , --- 


Ki+1 + Axk; + Bye ki-1 + Ci-2Ki-2 + ides 


, D; , E, , +++ , are also found to 


means of which the desired orthogonal rows and 


columns are generated thus become 


Kiss 


Ke+1 


uk; — ask; — Byki , (9) 


K;U — ayes Bi-1Ki-1 , (10) 
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where 
xuk; 8 K;- uk; 
es 9 ae se 
: kik; Ky-1k yy 


The pertinent scalars and zeros having been established, it is now convenient to identify 
the elements of the reduced matrix, outlined as follows: 


xwttky —xyttky nitty —_nytukg rity 











kk, Kk, eles Kk, Kk, 
Kuk, Kuk, Kuk, Kuk, K2UK 
Koke Koko Koko Koko Keke 
| | = | Kuk, Ksuk, xyuks xyuk, KjUK ny (11) 
kjk; - k3k3 k3ks3 k3ks3 k3ks3 K3ks 


Kuk, Kuk Kus Kuk, K UK 


Kyks tales Kyks Kak, Kyky 


KmUk, KmttKo Km, Kmuky Knltkm 


mm Ce iS KnKkm Km!Cm 

















The elements are easily recognized and replaced, yielding 


a B O DO «+ O 


- 0 1 a B -: O 
E ‘|= (12) 


Kk; 














0 0 0 O =e a, 


The vanishing of so many elements makes the expansion of the characteristic determi- 
nant extremely simple, so that a treatment of the reduced matrix equation by conven- 
tional iterative methods is unnecessary, and the step from the formation of the orthogonal 
row and column sequence to the solution for roots and modal columns may be taken 
without specifically examining the reduced matrix. The determinant will be expanded 
in stages, beginning with first order, which corresponds to the choice of one row and 
column, and successively increasing the order until the desired number of rows and 
columns has been reached. The resulting polynomials, individually equated to zero, will 
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have roots which will represent successive approximations to the roots of the original 
matrix. These polynomials, which may be written by inspection, follow: 


The determinant is 


JA — a —B, 0 tee 


| Tk —] A - A> — Bo — 7 
pia E | i = Dn(d). (13) 


The polynomials are 


D(A) =A- ay 


p2(X) = (A -— Q>)p,(A) — B, 
p3(A) = (A — as)p(A) — Bop,(A), (14) 
ps(A) = (A — as) p3(A) — Bsp.(A), 


Diail(A) = (A — ayai)pi(A) — Bipi-s(A). 


The solution for the modal column corresponding to a root of p,,(A) = 0 will be 


available after the coefficients, c; , have been obtained by direct solution of the reduced 
matrix equation as a set of simultaneous scalar equations. Since one coefficient may be 





arbitrary, let c, = 1. Then, 
A— a, p,(A) 
Co — a —_— = . 
B, B, 
(A — av)Co — 1 p(X) 
t — - — = = = = 
B 8,8 
(A = —_ (A) 
ae a Qs )C C2 oz P3\4) 
3 8, BoB; 
(A — ae; — Ci-; pir 
rf 4 — = = = —_— ~ ——_ 
. B BB» B, 


and we have the modal column through its initial definition, 


a (A) (X) m—1(A) 
k= Dek: =k, +?“ bk, + OY + ee + Pe. (18) 
B, B,B> B:B>o °** Bm 


It is important to note that, not only is this a finite series, but it is limited in length 
to a maximum number of terms, m = n, the order of the original matrix. The generation 
of the orthogonal function sequence stops at this point, since the original matrix on 
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which it is based is incapable of defining more than n linearly independent functions, 
and when this number has been reached the reduced matrix will yield precisely the 
roots and modal columns of the original matrix. This characteristic, more fully dis- 
cussed by Lanczos [2], is an important feature of the method of minimized iterations, 
since it indicates that this method, unlike the conventionai matrix iteration algorithm, 
will never require a continuation ad infinitum. 

4. Minimized iterations for the non-homogeneous equation. Given the equation 


(AI — u)k = q, (16) 


where \ and g are specified and a solution for k is desired, the reduced equation is formed 
exactly as for the homogeneous equation, replacing k and premultiplying by [x;] as 


(x7 — «aut: |). = | Jo. 
Kjh;j Kjhj 


The column on the right side of the equation is vastly simplified by letting k, = q 


(x7 _ Ea _ {1, 0, 0, ety 0}. (17) 


The method of minimized iterations is now employed to determine the same scalars, 
3, , as before, whereupon the solution for c is indicated as follows: 


a; and £; , 
Kuk; ee 
c=\|AIT —- ae {1, 0, 0, --+ O}. (18) 


The inverse of a matrix is defined as the quotient of the adjoint of the matrix by its 
determinant, and the adjoint is the transpose of a matrix whose elements are the co- 
factors of the corresponding elements in the matrix to be inverted. The determinant 
of the matrix is already available as p,,(A), so it is required to examine the adjoint. 
Since the adjoint is used to premultiply a column whose elements all vanish except the 
first, only the first row of the adjoint is needed. The elements of this row will be the 
cofactors of the elements of the first column in the determinant array of the reduced 
equation, (13). These cofactors can be written directly by inspection, using polynomials 
of the form already defined in Eq. (14) but in reverse order, beginning with the lower 
right corner of the array instead of the upper left. These reversed polynomials will be 


follows: 


whence 


designated as 
Di(A) =rA— ay ’ 
p2(A) = (A — Qm—1)Pi(A) = ’ 
(19) 
PsA) _ (A a Q@m-2)P2(A) saree Bn-2p,(A), 


Divi) = (A — Om-i)PiA) — Bui. 


Using these, the first cofactor is p,-,(A), the second is 8,p,,-2(A), the third is 8,8.),-3(A), 


and the coefficients, c; , may be written, 
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Dm— (A) 
¢, = Bau) 


Du(d) ’ 
a = Bs 
Bea oy 


Applying these to the definition of k, the solution is complete. 
- bas , , 
k= z ek; = p (r) (Dm—1(A)ky + B:Dm-2(A) ke + °°: + BBs *°* Burka). (20) 
i=l m\A, 


5. An outline of the Lanczos algorithm. While the preceding derivation may, appear 
to be complicated, practical application of minimized iterations to numerical calculations 
requires very simple repeated steps, which are outlined below and which may easily be 


translated into operations on automatic computing machinery. 
First stage: Given x, and k, . 
Form x,u, uk, , xk; 
Form «,uk;, 


. Kuk, 
Compute a, = ——— - 
Kk, 


p(A) =A —a, =0 
now yields a first approximation to the first root. 
Second stage: Form xk. = xu — a,k; , kz = uk, — ak, 


Form ott, uke , Keke , Kottke , ntke « 


Compute a. = ——, A; = 
Koke K, hy 


Kok Kuk, 
pA) = (A — a)p,(A) — 8, = 0 
gives a second approximation to the first root and a first approximation to the second. 


— Bk, . 


Third stage: Form x3 = ku — aok. — Bik, , ky = uke — ake 


Form x3, uks , ksks , Kglks , Kollks . 














1951] SOLUTION OF THE MATRIX EIGENVALUE PROBLEM 25 
uk; xouk 

Compute _ ae =. 

— ™ K3ks ‘ Kok» 


gives the third approximation to the first root, the second approximation to the second 
root, and a first approximation to the third root. The calculation is continued as far 
as necessary in continued stages. While the process of minimized iterations as described 
here can be expected eventually to provide the desired latent roots to any desired degree 
of accuracy with any initially chosen row and column, the rate of convergence can be 
made extremely rapid by first multiplying the initial row and column several times by 
the matrix, u, and then beginning the application of the method. The advantage of 
doing this is that the minimized iterations will then begin with functions whose normal 
mode components will be “ordered” in the sense that the dominant mode will represent 
a large part of the initial functions, the second mode will be somewhat smaller but next 
in importance, the third mode will follow, and so on. To a rough degree of approximation, 
the orthogonal functions thus generated will then in turn approximate the normal 
modes in the same order, and the roots of the successive polynomials will also appear 
in this order. Lanczos [2] presents an impressive numerical example which demonstrates 
extremely powerful convergence. 

It is convenient, particularly in the use of punched card machines, to form certain 
of the required products simultaneously and to supplement the required operations by 
check calculations which verify the orthogonality of the generated rows and columns. 
As an example, consider the third stage. The product, x3u, is indicated first and the 
matrix is bordered by several columns as follows: 


K3[U, ky ’ kp ’ ks] _ [k3u, 0, 0, Kk]. 


The two zeroes serve as a numerical check. Likewise, the product, uk, , is also formed 
in conjunction with others by bordering u with several rows, as follows: 


[u, K,U, Kou, kzulk,; = {uk , 0, xouks , xguks}. 


The zero here serves also as a check, and all the scalar products required for the com- 
putation of a; and 8, are available from these operations. In numerical calculations, 
when the zeroes may occur as small finite numbers due to rounding errors, a refinement 
is possible by using these numbers to determine the constants, y; , 5, , €; , ete. of Eq. (7) 
and (8), which may be employed as indicated by these equations for the generation of 
further orthogonal functions. 

6. Minimized iterations in the Galerkin method. It will now be shown that a mini- 
mized iteration technique, similar in form but slightly different in detail, is also applicable 
to the solution of the eigenvalue problem by the Galerkin method, and that this method 
may offer certain computational advantages. Beginning with the reduced equation (5) 


i 


in Galerkin form, but with kik; = 0 for 7 ¥ J, 


Kjuks |). 
(7 a Es } — (21) 


the generation of suitable columns, k, , will be detailed. Using the same technique 


proposed by Lanczos, let 
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heey = uk; — ak; — Byaykgny — Yi-aki-2 — 6:-3hi~2 °° 
and solve for the constants required on the basis that k¢k, 0 (i ot j). 
h Rix = kiuk, — a, kik; _ 0. 
kikise = kyukisr — Bikik; = 0, 
kik; Keful - y kik; = 0 
whence 
a a kilo kiuk +2 
ke} P p; = a ’ gf I'he . ete. 


However, in the Galerkin case none of these constants vanish. Nevertheless, 


| A wk = k: h — 2 
and 
ki. uk; = 0 if yo t. 


Thus, the reduced matrix becomes 


Ea _ 0 l As Bs =o 1 
ik, || . iP 


which is of the same form as Eq. (11), (compare with Eq. (12)), and the 
expansions of the determinant by steps are, beginning with the definition, 


| Am —B, sa 


| 
—1 
kiuk, | 
a7 a a od 


A — ae —B 
= Pmid); 





p(t) =A—a,, 
P(A) = (A — aa)pi(A) — B , 

p3(A) = (A — as)p(A) — Bopi(A) — 1 » 

pA) = (A — ag)ps(A) — Bspo(A) — yopi(A) — 6, 


etc. 
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(22) 


polynomial 
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The roots of these polynomials, as in the Lanczos method, yield successive approxima- 
tions to the latent roots of the original matrix. 

The solution for the coefficients, c; , required in order to define the modal columns, 
must here be accomplished in reverse order, due to the triangular form of the reduced 
matrix. These coefficients will be expressed in terms of a set of reversed polynomials, 


defined as follows: 


Ca-1 = A — a, = D(A), 
Cc o> (A — Am )Cm—1 a | — (A — Qm—1)pi(A) ual Bini = p.(a), 


Ca-3 = (A _ Qm—2)Cm—2 os B.. Pan Ym-2 ’ 


= (A — an 2)p2(A) v3 Bm—2p(A) = To~— * ps(A), 


etc. 


Hence, the modal column is written for \ corresponding to a root of p,,(A) = 0, in terms 
of the reversed polynomials defined above. 


k = Decks = Pmald)ky + Pu-a(ho + +++ + Phen + ken - (26) 
=1 
It is interesting to compare Eqs. (21)-(26) with their corresponding relationships in the 
Lanczos method, Eqs. (6)-(15). 

7. The Galerkin treatment of the non-homogeneous equation. In a manner similar to 
that applied by Lanczos and described in Sec. 4 for the solution of the non-homogeneous 
equation by minimized iterations, the method of minimized iterations may be extended 
to the solution of the reduced equation obtained by the Galerkin method. The necessary 
scalars are determined as described in Sec. 6, thus permitting the determination of the 
two sequences of polynomials, p;(A) and p,(A). The next step is to solve the reduced 
equation for the coefficients, c; . To accomplish this by determining the cofactors re- 
quired from the adjoint of the matrix to be inverted is not convenient, in view of the 
triangular matrix form, but the solution can be easily obtained by first solving for the 
coefficients in terms of the last, c,,, working upward in the array until the use of m — 1 
equations provide c, through c,,-, in terms of ¢,, , and then using the first equation to 
evaluate c,, . This process follows: 


Cm-1 = (A — Gn)Cm = DilAlCm , 

Cma2 = (XN — Gm-1)Cm-1 — Bm-iem = P2lA)Cm ; 

Cu-3 = (A — Gn-2)Cm-2 — Be-2lm-1 — Yun-2Cm = DalAlCm » 
etc. 


[(A _ 1) Pm—1(A) _ BiPm-2(A) i V1Pm-s(A) “* * Be = 1. 
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The coefficients having thus been determined, the solution follows directly. 


as =; ie Du (Nk, + D, AAN)ko + e+ + D(Nkn- + ke 27) 
“1 (X — @)Pm—-1(A) — BiPm—2(A) — V¥iDm-a(A) *°° 


(27) does not closely re- 


In view of the means by which the c, were determined, tq. 
semble Eq. (20). 

8. An outline of the minimized iteration algorithm in the Galerkin method. Again it is 
desirable to present an outline of the actual steps in numerical calculation by minimized 
iterations, from which a logical sequence of operations may be established for purposes 


of automatic calculation. 
First stage: Givenk, . 
Form uk, ; kik, . 
Form ktuk; 


Cc kiuk, 
ompute a, = zn 
kk, 


pi(A) =A—aQ,. 
Second stage: Form k, = uk, — ak, . 


Form uk. , kik. , kzuk. , kiuke «3 


& t kzuk, 3 kiuk, 
ompute a, = —  §=—-- 
‘i kok, Kak 


plA) = (A — ae)pi(A) — Bi 
Third stage: Form k, = uk, — asks — Bk, 
Form uke . kik; , k ‘uke : kiuk, . kiuks ° 


keluk kiuk, kiuks 


Compute =— fo = <= 3 = 
ee ee 


Dp (A) = (A -—a pArA) — B.p,(A) a 


As in the Lanezos procedure described in Sec. 5, it is advantageous to precede this proc- 
ess with several premultiplications of an arbitrary column in order to control the 
relative magnitudes of the various modal components of the column, k, , so that the 
roots will be obtained from the successive polynomials in order of relative magnitudes. 
As indicated for the Lanczos algorithm, it is convenient in punched card procedures to 


form certain of the required products simultaneously and to include check calculations. 








; 
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To this end the matrix, u, is bordered by a sequence of rows, as exemplified below for 
the third stage. 


[u, 4 4 ’ k3lks - {uks ’ 0, 0, k3ks} ° 


The two zeroes serve as a check on the orthogonality of the generated columns. In 


« 


addition, the following operations are necessary: 


[es , ke , Keifuk, = (heluck, , keluk, , kiuks}. 
whence the constants, a, , 6. , and y, , may be calculated. 

9. General comments on minimized iterations. It should be noted that the foregoing 
presentation and extension of the method originally proposed by Lanczos has leaned 
toward the pictorial approach to a means of mathematical analysis rather than con- 
cerning itself with important considerations of mathematical rigor. The special cases 
of equal roots and of unfortunate choice of arbitrary row or column, where a deficiency 
in modal components might lead to peculiarities in the results, have been deliberately 
ignored, first, because these cases are adequately treated in Lanczos’ paper [2], and 
second, because they would detract from the object of presenting a practical procedure 
and of emphasizing the computational aspects which are of interest to the applied 
mathematician and to the engineer. Furthermore, in the interest of providing a simple 
derivation which would admit of obvious extension to the Galerkin type equation, the 
original intent of the name, minimized iterations, was passed over lightly, deriving the 
basic scalars through strictly algebraic operations without consideration of their more 
fundamental purpose of providing the maximum utility with the minimum number of 
row and column functions. The advantages of orthogonal functions in computation 
provided an equally direct reason for the choice of these scalars. 

The method of minimized iterations, in either the Lanczos or the Galerkin varia- 
tions, is recommended as an alternative to conventional procedures in the numerical 
solution of matrix problems on the basis of potentially large time-saving in computation. 
Its advantages are most evident when the determination of a number of eigenvalues, or 
latent roots, is of prime interest, and the determination of modal columns or the solution 
of non-homogeneous equations can also be greatly expedited, subject to an efficient 
organization of problems for whatever type of computing machinery may be available. 

The Galerkin variation appears to be of greatest advantage in the eigenvalue prob- 
lem, requiring substantially only half the number of matrix-column products involved 
in the Lanczos procedure. However, the computation of modal columns by this method 
is somewhat more cumbersome, particularly if successive approximations to a modal 
column are to be studied, in view of the requirement for reversed polynomials. In the 
non-homogeneous equation, on the other hand, this disadvantage applies to both ver- 
sions, the only significant difference being that the Galerkin polynomials are longer. 
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SOURCE-SUPERPOSITION METHOD OF SOLUTION OF 
A PERIODICALLY OSCILLATING WING AT SUPERSONIC SPEEDS* 


BY 
H. J. STEWART anv TING-YI LI 
California Institute of Technology 


Introduction and summary. In a recent paper Evvard (Ref. 1) discussed the linearized 
theory of the non-steady motion of three dimensional wings by methods which he had 
previously developed for the treatment of the corresponding steady flow problems 

tefs. 2 and 3). Evvard represented the wing by a distribution of sources, and the 
important result of his steady state theory concerned the=letermination of the flow in 
region influenced by a subsonic leading edge or wing tip. He showed that the influence 


of the flow around this subsonic edge of a flat lifting wing on the velocity potential at 





x=X,(y) 
4 “P 0 _ 
; — x=%p( y) 
- = 
e pint 





x 


Fia. 1. Geometry of Wing. 


a point within the region of influence of this edge is exactly equal to and of opposite 
sign to the contribution to the potential from the sources distributed over a simply 
determined region of the wing. In his paper on non-steady motion, he was able by similar 
methods to determine an explicit formula for the velocity potential; however he could 
not express the results in a similar, ‘equivalent area’’, form. 

The present paper is concerned with the same problem of the non-steady lift of 
finite wings at supersonic speeds, particularly in regions which are influenced by subsonic 
leading edges or wing tips. It is shown that the simple “equivalent area” theorem de- 
veloped by Evvard for the steady state case is also valid for oscillating wings. The 
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theorem is not extended to arbitrary non-steady motions, and an example where the 
theorem in this simple form is apparently not valid is demonstrated. 

The basic differential equation and boundary conditions. Consider a wing in a steady 
supersonic flow of velocity U and Mach number M in the direction of the x-axis. Then 
the velocity potential ¢ which governs any small, possibly non-steady, disturbance pro- 
duced by the wing satisfies the linearized differential equation 


BF 


a’ at a odxoat 


a 


1 dg | QU dy ( i 1) de _ de , He (1) 
On” Oy” d2”’ 


where (x, y, z) are cartesian coordinates, ¢t is the time variable and a is the speed of sound 
in the undisturbed flow so that U = Ma. The wing is assumed to be near the plane 
z= 0. 

With the approximation of the linear equation, it is permissible to replace the 
boundary conditions at the point (z, y, z) on the actual wing surface by the same boundary 
conditions applied on the plane z = 0 at the point (x, y, 0). In order to express the 
boundary conditions it is necessary in general to divide the wing surface into two different 
types of regions (see Fig. 1). The origin of the coordinate system is taken at the point 
0 where the Mach line Ov is tangent to the leading edge. The leading edge is thus divided 
into two segments, the segment 0A which is defined by x = 2,(y) or y = y,(x) and is 
a supersonic leading edge and the segment 0B which is defined by x = x,(y) or y = y2(2) 
and is a subsonic leading edge. As a matter of convenience it is assumed that the trailing 
edge, x = 2;(y), is a supersonic trailing edge where the Kutta condition need not apply. 
The Mach line Ou then divides the wing into two regions. Region I, which is bounded 
by x = 2,(y), x = 23(y) and Ou, may be referred to as a purely supersonic region. Region 
II, which is bounded by x = 2.(y), x = 2x3(y) and Ou may be referred to as a mixed 
supersonic region (Ref. 4). 

At any point on the surface of the wing the flow must be tangential to the surface 
at any instant. This boundary condition, linearized, and applied to an oscillating condi- 


tion is 
Np 7 , ; ‘ . 
= = w(x, y, + 0) exp (ivt) = UA,(x, y, + 0) exp (ir), (2) 
dz 
Ip , . . 
= = w,(x, y, — 0) exp (vt) = —UA,(z, y, — 0) exp (i?), (2a) 


where, except for the time factor, w(z, y, z) is the z component of the velocity and A 
is the effective slope of the streamline and » is the frequency of oscillation. The subscript 
T refers to the top of the wing and the subscript B refers to the bottom of the wing. In 
general w, and w,; (or Ay and A;) need not be related. A sign convention for A; and Ag, 


adopted in Ref. 1 is also used here and is shown in Fig. 2. 


From the definition of a purely supersonic region, there can be no disturbance in the 
flow ahead of the line x = 2,(y). For any point P in Region I the velocity w is thus 


known at every point on the plane z = 0 in the forward Mach cone from the point P. 


On the wing w is given by Eq. 2 or Eq. 2a, ahead of the wing w = 0. 
tegion II conditions are more complex. As before, the velocity w 


For a point Q in 
is given for that portion of the plane z = 0 in the forward Mach cone from Q which is 
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covered by the wing by Eq. 2 or 2a. Also w = 0 and g = 0 ahead of the line segments 
x = 2,(y) and Ov. Since x = 2,(y) is a subsonic leading edge, there is, in general, an in- 
teraction between the upper and lower surfaces which produces an upwash across the 
plane z = 0 in Region III which is bounded by « = 2,(y) and Ov. This upwash cannot, 
in general, be specified in advance. For this region the pressure must be continuous 


across the plane z = 0 so the linearized boundary condition for this region is thus 
Oger , 77 9er _ In | 77, We , 
a ee oe tO @) 


The boundary conditions on the plane z = 0 for a point Q in Region II are thus of a 
mixed type, involving w over the wing, pressure over Region III and no disturbance 


ahead of the lines Ov and x = 2,(y). 








Fic. 2. Sign convention of A’s. 


Elementary oscillating source potential. Elementary solutions of Eq. 1 which can be 
superimposed to obtain more complex solutions can easily be obtained by the method 
f separation of variables. For this purpose it is convenient to introduce the following 


ol 


coordinate transformation: 


r=(— By +2)”, 


| B(y? + ay 
pliner 


\ 2 A 
w = tan’ (z/y), = sal of v2), 


(4) 


where 8° = M*® — 1. 


These space variables were found useful in the treatment of steady conical flows (Ref. 5). 
The time transformation is similar to a combined Lorentz and Galilean transformation 
and has been used by Miles (Ref. 6). In these coordinates Eq. 1 is 


Oe »0¢ 0g 0 | ° ge | ] Fe) 
, =r 2r - —!i(Q1—- —— ess 
" Or or” sas or + Ou v) Ou ¥ 1 — yp dw (5) 


This is identical with the form of the classical wave equation in spherical coordinates. 
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The solutions of Eq. 5, obtained by the method of separation of variables are 


( cos ms (PZ(u)) (r : = 1 a(lr)) 


e= 2) Arm | 
se sin a lg u)) B I anrhl) j 


where /, m and n are the separation parameters. P*(u) and Q?"(u) are Associated Legendre 


exp (+7l7), (6) 


functions and J «;,+:,/2) (lr) is a Bessel function of order +(n + 1/2). 
For m = n = 0, a simple solution of Eq. 6 is 
g, = Ar’ J_,,2(lr) exp (ilr). (7) 


If + is replaced by the physical time variable from Eq. 4, the Bessel function is written 


in the trigonometric form and 1 is eliminated by the relation v = /6°a, Eq. 7 becomes 


Rae cos a) exp E {-— Par | (8) 


where A, is a new arbitrary constant. Equation 8 may be considered as defining a super- 
sonic oscillating source. This basic solution has been used in this form by Miles (Ref. 7). 
If the basic solutions used by Garrick and Rubinow (Ref. 4) or by Evvard (Ref. 1) are 
applied to oscillating problems, they can be reduced to this same form. It may be noted 
that for y = 0 Eq. 8 reduces to the usual steady state source potential. The complete 
velocity potential field for an oscillating source is defined by Eq. 8 in the downstream 
Mach cone and as zero outside the Mach cone. 

Velocity potential of an oscillating wing. For a point in the purely supersonic region 
the velocity potential due to the wing can readily be obtained by replacing the wing 


Zz 





x,Y,Z) 





(x-ez,y,0) 


Fic. 3. Singularities or sources in the z, y plane, that affect conditions at (2, y, z) at instant t. 


by a distribution of sources over the wing surface. If the region of dependence of a given 
point includes only that portion of the wing which is purely supersonic, the velocity 


potential for z > 0 due to the source distribution is thus 


i f U(x — &) \ vr, \ dé dy 
, te = E, - S a . Ss ] F 
&(z, y, z, t) = [| A7(é, n) exp jw} t — Ba? f cos (4) . (9) 


Je ( 
where 
2 2/ ‘s 92,241/2 
r= {[(#—-— —-BYy-—n) —- Be”. (10) 
Here A-,(é, 7) is the source strength per unit area at the coordinate (&, ») on the wing 
surface. The region of dependence, which determines the region of integration on the 
wing surface, is bounded on the downstream side by the line 7, = 0 (Fig. 3). 
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If the source strength A,(é, 7) can be chosen so the boundary conditions for the 
purely supersonic region are satisfied, Eq. 9 is the proper solution. In order to do this 
it is convenient to replace the integration variable y by 


B(y — n) = —r, sin 8, (11) 
where 
Tr. = [(x — &)? — B’e"]'”. (12) 


With this notation Eq. 9 becomes 


r—Bz ° 7 
P(x, y, z, t) = 5 OP (ivt) £ exp | -80, (x -—é | dé 


(13) 
(~ ( T. . ) ( vr, 46 
a A7r\é, yt ame cos Ba cos @ : 
where £, is the least value of — on the leading edge. Since 
Or, _ _ Bz 
de r? (14) 
“* = —2A,(x — Bz, y) exp (1) 
dz 
— ge exp (ir) [exp | 20 — | hae (15) 
Bz exp (iv I. xX} Ba x—¢ z 


eek _— UIs 
Bs Je 2 290s |— cos E 
* | ars E (e, y+ 3 sin 0) COs (#: co 6) | dé 


If the function A ;(£, 7) is continuous in the neighborhood of the point (x, y), the magni- 


rary 


tude of the double integral in Eq. 15 is finite for sufficiently small values of z; so 


lim (22) = —A,(z, y) exp (iv). (16) 


Oz 


By comparison of Eqs. 16 and 2, it is seen that 
I U 
A,(z, y) = —-— w(x, y, + 0) = —— A(z, y, + 0). (17) 
T T 


For a point below the wing, z < 0, a similar analysis shows that 


U : 
A,(z, y) = J w,(z, y, — 0) = —— A;(z, y, — 0) (17a) 
T 1 


The required source strength, A(z, y) in the plane z = 0 is thus completely determined 
for any point in the purely supersonic region. On the wing A(z, y) is given by Eq. 17 or 
Eq. 17a. Ahead of the wing the disturbance (0@/0z),-) is zero; so A(x, y) = 0 in this 
region. With these values of A(x, y), Eq. 9 defines the velocity potential and thus the 
velocity components and the pressure on the wing. This analysis was given in a similar 
form by Miles in Ref. 6 (some errors in his presentation were corrected in Ref. 7). 
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The results of this analysis may be summarized in the following two theorems: 


Theorem 1. The strength of the source at any point at any instant on the surface of 


an oscillating wing is linearly dependent on the downwash at that point and at that 
instant and is independent of the downwash of the neighboring points. 
Theorem 2. The velocity potential at instant ¢, at a point P in the purely supersonic 


region of the surface of a three-dimensional oscillating wing (Fig. 1) may be computed by 


; U : ” ~ 7 
Hz, y, + VU, t) = —— exp wt) | exp | —tv s3-3 (a — &) | dé 
T Je, pa 
(9a) 
nyt (2 gs {Ar(E, n) 2 2 271/2 
: at j sails ( cos {(v/Ba)[(x — &)° — Bly — n) |'”"} 
x | ) | _— (xz — $2 ia e*(y — 9) |" ay, 
ey Th ae ‘ins 
; P top ‘ . ‘ 
where z = +0 refers to the surface of the wing. 
bottom 


A mixed supersonic region may be converted into a “‘psuedo-purely supersonic region”’ 
by Evvard’s procedure of inserting a diaphragm into Region III of Fig. 1 which is an 
extension of the wing having the following properties: 


a) It does not change the flow over the wing 


b) It sustains no lift 
With this supposition, the top and bottom surfaces of the wing may again be considered 
to be independent so that Eq. 9a applies; however the diaphragm slope is in general an 


unknown function which must be determined. 
A part of the wing in Fig. 1 is shown enlarged in Fig. 4. In order to compute the 





U 
Ve ty) 
\ X= X,(y) 
I 0 y 
(Xp, Yp) 
TSX (uo.te) 
(uw) ._H Vv 
Q \_ Ve Ve (u) 
u (x,y) X= Xo(y) 
| (Uw, Vw) 





x 
Fig. 4. A portion of the wing in Figure 1 


velocity potential at the point Q, it is convenient to first consider a point D located on 
the trace of the upstream facing Mach cone from Q, in the diaphragm plane. Let the 
unknown downwash and the effective slope of the streamline on the top surface of the 
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1951 
diaphragm be wp,(xp , yp) and Ap,(Xp , Yn) respectively. Then, by Eqs. 9 and 17, it is 
found that 


l 7 
On (Xp, Yn, + O, t) = —— exp (a) 
T 
- — 12 ———— = 2z1/2 5 :: 
x [| este: = 1 6 a)((a D ~§) : B (yp — 1) | “ — dé dy 
J. tp — &) — B(yo — 0) | exp |[W(U/B'a)(an — &)] 
(18) 
_ : exp (ivt 
' cos {(v/B’a)[(x» — &)” — B(yn — n)|'”"} 
x [ u (€ 7 ~ a TN © - 4 y maa \2 72 ae sar fae se d d 
] ‘ (tp — £)? — B(yp — n)|'” exp [(U/B'a)(ap — &)] £ dn, 


where Sy is the region of the wing and S, is the region of the diaphragm included in 
the upstream facing Mach cone from D(zxp , yo , + 0), at instant ¢. The regions of inte- 
gration Sw and S, are most easily expressed in terms of the oblique u, v coordinates 


defined as follows: 


M : 8B 
oe a 
u 28 (é Bn), E M (v + u), 
or (19) 
M 1 
v= 28 ( + Bn), » eae M (v — u), 
With these coordinate transformations, the point (xp , yp) is transformed into (up , Up), 
where 
M B 
Un = 95 (tp — BYp) in = 7 (vp + Un) 
Oo (20) 
iy i es 1 hit. ~ tJ 
o = 28 “ D Yb Yd M D D 


The surface integral @p, in Eq. 18 will now be integrated in the u, v plane. Then, 


Eq. 18 becomes 








l . x du 
os .. +008 -—e _ 
Porun 1%, +9, sm “P (wt) Jo (Up — u)'” exp [(tv/Ba)(up — u)] 
- = Wr(u, 2) cos {(2» ‘MBa)|(up - ud(vp — v)]'"} ip 
Jester (vp — v)'” exp [(iv/Ba)(vp — v)] 


(21) 


an —s * 7 du 
5 Y ies ) Jo (up — u)'” exp [(tv/Ba)(up — u)] 





et 


x [~ Wp,(u, v) cos {(2p ‘MBa)|(un — u)n — v)) 
(vp — v)'” exp [(iv/Ba)(vp — v)] 


“valu 


dv. 








where w,(u, v) is the downwash on the top surface of the wing, wp,(u, v) is the down- 
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wash on the top surface of the diaphragm, the area bounded by 0 < u < up and v,(u) < 
v < v.(u) is Sy , and the area bounded by 0 < u < up and v,(u) < v < Up is Sp (Fig. 4). 

Similarly, for the corresponding point D(up , vp , —9, t) on the bottom surface of 
the diaphragm, it is seen that 


®,,( 0, ¢ : (ivt [~ _ 
(up , Vp, — 0, 2) = > exp (at) EMS 1. ORE TCR Oa? NE 
ee ee rM “Jo (tp — uu)” exp |(tv/Ba)(up — u) 
w,(u, v) cos {(2vy/MBa)[(up — u)(vp — v)]'/*} 
x | . ta |. ee , tr amos Ty — dv 
ee (vp — v) ~~ exp [(tv/Ba)(vp — v)| 
(22) 
1 puD du 
+. exp (art Re Fe FIG ET ES 
rM , (up — u)’ exp {(2v/Ba)(up — u)] 
x [ > wWp,(u, v) eos {(2v/MBa)[(up — w(vp — v)|'"*} J 
: ae — nae meme 
| vp — v)'”” exp [(tv/Ba)(vp — v)] 


Off the wing, the downwash must be continuous. In terms of the effective slopes of 
the stream lines, this condition is, with the sign convention of Fig. 2, 
Ap,(u, v) = —Ap,(u, v) = Apu, v). (23) 
From Eq. 3 it is found that in Region III 
},,(z, y, + 0, t) = Bp,(z, y, — 0, 2) + F(a — Ut, y), (24) 
where F is an integration function. The foremost Mach cone (Fig. 4) from the origin, 
0, represents a line of infinitesimal disturbance along which F(x — Ut, y) can be set 
equal to zero at all times. F remains zero along y = constant lines for values of x not in- 
tercepted by the wing (Ref. 1). Therefore, in Eq. 24, F may be put to zero. Behind a 
trailing edge, F may be different from zero and the theory must be modified. Then, from 
Eqs. 21, 22, 23 and 24 (with F = 0), it is seen that 


] [ du 
2 “ Un = U ‘, exp [(av Ba)(Up — u) | 
A [ - A,(u, v) — Aru, v) | cos {(2v/MBa)|(up — wlvp — v)}'”7} ] 
Pe — — = . i > 7 . aa on av 
oe Wn = U) exp [(ay ba)(vp — v) | 
(25) 
[ du 
ma lp — u)'” exp [(iv/Ba)(up — u)] 
eat ? Ap(u, v) cos }(2y MBa)[(up — u (vp — v)]'”*} 
x 7 ay dv. 
(vp — v)~ exp [(tv/Ba)(vp — v)| 
When v = 0, this reduces to 
SP i du nos Ap(u, v) — Aru, v) 
—— 72 _ 72 dv 
23006 (Up — (vp — v) 
(26) 
a uD du avD Ap(u, v) 
= re 172 / 17a dv. 
0 (Up — U)* Jou) Ur — Y) 
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Inasmuch as the limits of integration of the u-integrals are the same for all values of 
up and the integrals are “faltung’”’ integrals, the two integrals with respect to v may 
be equated along lines of constant u that extends across the wing and the diaphragm 
(Fig. 4). Therefore, 


f° Ante, dy =f" Avle,) — Aele 9 a (27) 


(vp — v)'” Fiutus 2vp — v)'” 

This is the fundamental result of Ref. 2 and also is the basic equation of Ref. 3. The 
above argument is valid because the terms containing up do not appear in the v-integrals, 
and hence Eq. 27 is true for all aie on the line v = pvp. 

The parallel treatment of Eq. 25 would be possible if the terms containing (up — u) 
(vp — v) can be separated as in Eq. 26, under the integral signs. The present treatment 
represents a first attempt towards this end. The isolation of terms containing (up — u) 
from terms containing (vp — v) such that the v-integral is free of the (up — u) factor, 
may be accomplished by the following procedures. 

The term {(vp — v) (up — u)}'” vanishes at (up , vp), therefore Eq. 25 actually 


should be 


puD—e 


lim | , 
(ty 


e-0 “0 ) 


du 
— u)'” exp [(iv/Ba)(up — u)] 





x [ um [Antu, oa Ay re _v)] cos {(2»/MBa)[(up — wn — “ay 
2» — v)'* exp [(ix/Ba)(0» — »)] 


(25a) 


lit = ae Ss ee du 
‘ . ‘ (up — u)'”? exp \(wv/Ba)(up — - u)] 


dv. 





fs re" A,(u, v) cos {(2y/ MBa)[(un — wp — v)}'/7} 
i (vp — v)'” exp [(tv/Ba)(vp — v)] 


The nature of the functions Ag , Ar and Ap must be such as to insure the existence of 
the improper integrals. Thus, except for the singularity (wp , vp); in the finite integration 
regions, the integrands are defined and bounded everywhere. Now, the circular functions 
are defined by power series; in particular, the power series expansion of the cosine func- 
tion 1s 


COS z « ie sae . (28) 


The series (28) has the following properties (Ref. 8): 


(1) It converges absolutely for all values of z (real and complex), 
) It converges uniformly in any bounded domain of values of z, and consequently, 
It is a continuous function of z for all values of z 


(2 
(3) 

Because of the uniform continuity, the cosine function in Eq. 25a may be expanded 
in an infinite series and the orders of integration and summation may be i verted. 


Thus 
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) 


x Ii l = =a x " TAp(u, v) — Ar(u, v) (vp — v)"'” j 
im - A ___——— = Da fm A - dy 
paren & exp [(7v/Ba)(up — U)] Jose 2 exp [(7v/Ba)(vp — v)] 
(25b) 

yd —)"(vy/Ba)?"(1/M)""2'”” 

 & n'iT(n + 1/2) 
a —€ n-1/2 at —€1 \ n—-1/2 

; ln — U) du Ap(u, v)(vp — v) 
x lim | . a | pcre? moet: a, 

ere exp |(zv/Ba)(Up — U)] «,-0 .) exp [(tv/Ba)\Up — v 

I 1 2 rt 


With the conviction that the improper integrals under question exist, the “lim” signs 
may be left out. 

In Eq. 25b, unlike in Eq. 25a, the v-integrals do not contain up terms, and the problem 
has been reduced to one analogous to that of Eq. 26. Now, it may be pointed out that 
since Eq. 25a is derived by equating the velocity potential on the top and bottom sur- 
faces of the diaphragm in Region III (Fig. 4), the two sides of Eq. 25b may be con- 
veniently considered as power series in (1/1) of a potential function %, satisfying the 
original linear differential equation, Eq. 1; consequently corresponding terms may be 
equated. 

Therefore, for constant value of vp , with n being any positive integer, 


erp ant 4 J f / a . : _ ,\" 4 
An(u. v)(vp — v) | Ap(u, v) Ar(u, v)J&n_ - v)" wm (29) 
Jescu) exp [(2v/Ba)(vp — v)] 2 ee 2 exp |(2v/Ba)(vp — v)] 


In this system of simultaneous integral equations A,(u, v) and A;(u, v) are known 
while Ap(u, v) is unknown. Consider, say, (NV + 1) integral equations corresponding to 
n = 0,1, 2, --- N. (Of course, in the limit, VN —@). In order that these (NV + 1) simul- 
taneous equations may determine one unknown Ap , it is necessary that the (NV + 1) 
equations are not mutually independent, that is, the (V + 1) equations are reducible 
to one equation. In fact, this is true for the given system. For instance, when n = 1, 
it is obtained from Eq. 29 that 


r°e> ADU, v)(vp — v)'”” dv pra“) TA B(u, v) — Ar(u, v) (vp — v)'”? di al 
Sef oe | p—7_“". (30) 
2 exp [(zv/Ba)(vp — v)] 


exp (ap ba)\Up — UV) 


Carry out a differentiation of Eq. 30 with respect to vp . The result of this differentiation 


plus (iv/Ba) times Eq. 30 yields 


[~ Ap(u, v)(vp — v)'”? dv [ 7 TAn(u,v) — Ar(u, v) (vp — v)7'”* dh 31) 
— = _ (3 
exp [(iv/Ba)(vp — v)] 2 exp [(zv/Ba)(vp — v)] 
which is Eq. 29 for n = 0. Therefore, when Ap satisfies Eq. 30, it also satisfies Eq. 31. 


This argument can be carried on, by induction, to include the case for every n. Therefore, 
the system given by Eq. 29 is consistent and determines an unique function Ap . 
For the determination of the contribution of the diaphragm on the velocity potential 


at a point Q + (0) on the top or bottom surface of the wing, it is not necessary 
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to solve the integral equation, Eq. 29, explicitly. Let this contribution be called 
Dy,(Uw , Uw, + O, t) (see Fig. 4). Then, 


r 


l du 
Py (lw ,lw, +90, = F- t [ = BoE 
- nil qo _ Jo (uw = — 4)” * exp [(iv/Ba)(uw — u)] 








” Ap(u, v) cos | (2y ‘M Ba) [(uw om uw ak v) 7} dv 
(vy — v)'”” exp [(tr/Ba)(vw — v)] 
(32) 


4 U (—)"v ‘Ba)*"( “i/ ‘M)"*2'”” 
7 gre exp (tt) >! nIT( (n+ 1/2) 


0 


E on eile [" Ap(u, v)(vw — i dy 
exp [(tv/Ba)(vy — v)] ’ 


exp \(av/Ba)(uy — u) | di 


where u’ is the u-coordinate of the intersection point of the curves: v = v.(u) and v 
By comparing 


with — 


; Ap(u, v)(v_y — v)"'” ; pr’? Apu, v)(vp — v)"" ses 
| : dv, 
exp IG Ba)(vp — v)] 


exp [(7v/Ba)(v, — v)] 


it is seen that they are identical if every vp in the latter is replaced by vw . But the value 


of vp» along the v = constant line passing through the point (uw , vw , + 0) is vw (Fig. 4). 
Hence v, may be replaced by vw in Eq. 29 and Eq. 32 becomes 
U (—)"(v/Ba)?"(1/M)""2'”* 


Py (Uy Vy, £0, OY = FF exp (vl) orm — 
ua = nil(n + 1/2) 


n=0 


[ "(up — uy" '”? du [ mm [Ap(u, v) — Aru, v) Ow — eo yp 


exp [(2v/Ba)(uw — u)] Josao 2 exp [(tv/Ba)(vw — v)] om 
(33) 
= + e exp (iv?) ia } . 
aM J, (uw — ue exp [(iv/Ba)(uw — u)] 


i= [ *™ TAn(u, v) — Ar(u, v)] cos {(2v/MBa)[(uw — Wow — - v)]'” tuk 
Quy — v)'”” exp [(tv/Ba)(vw — v)] : 


In Eq. 33 an important theorem is established. The theorem may be stated as follows; 


Theorem 3. In the computation of the velocity potential at an instant ¢ at a point 
Q in the mixed supersonic region of an oscillating wing at supersonic speed, the contri- 
bution of the diaphragm may be evaluated by Eq. 33. In other words, the contribution 
of the diaphragm can be evaluated by an equivalent integration over a portion of the 
wing surface. Now, the velocity potential © at point Q on the top wing surface at instant 
t may be computed. It is (Fig. 4) 
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Duy ,uw, + O, 2) 


U Lig ~ du 
_—- exp (2rt) Ee 7 ae $$ —____—— 
7M Jy (Ug — u)'” exp [(tv/Ba)(uw — u)] 


Il 


*™ TAp(u, v) — Ar(u, v)] cos {(2»/MBa)[(uw — uw — v)]'”7} 
x . ee aa Seng NET EI er 


73 dv 
Vw — v)’~ exp [(tv/Ba)(vw — v)] 
U “1 du 
— — exp (ir?) ——— [Ne TS oy 
aM (uw — u)'’” exp [(tv/Ba)(uw — u)] = 
(34) 
oi f° ” A,(u, v) cos {(2v/MBa)((uw — u(ow — v)]'”7} 1 
: —— ———_—_—— — dy 
G (Vw — v)* exp [(vv/Ba)\Uw — v) | 
U ee r" ti du 
_ exp (ivf) ’ er ; : are a 
mM °™ (uw — u)'”” exp [(tv/Ba)(uw — u)] 
ce A,(u, v) cos {(2v/MBa)[(uw — uw — v)]'7} ] 
cae. we _ —— dv. 
(vw — v)’”” exp [(v/Ba)(vw — v)] 


In Eq. 34 the first surface integral represents the contribution from the diaphragm, 
while the last two surface integrals are the contribution from the top surface of the wing. 
By combining the first and second surface integrals, it is seen that 


Guy .vw, +0, 8) 


U ; ig du 
=-- exp (ivt) aaa See eae Pe -<_ 
rM (uw — u)'’” exp [(iv/Ba)(uw — u)] 
x [ y [A,p(u, v) + Ar(u, v)] cos {(2v0/MBa)[(uw — Wow — v)}'7} j 
4 2 - See ee a i72 ; ————— ee ae iia. 
oD as 2(v 4 — 7 exp [(av pa\vy — v) | (35) 
U : ‘aie du 
_— exp (2vl | iE FER — 
7M Jue (uw — u)’” exp [(tv/Ba)(uw — u)] 


v 2p / 1/2 
x | A u, 2 cos (Qy MBa)|(uy = U)(Vy ~~ v) | —- dy. 


1/2 (2 / 
Fa, (vw — v)’~ exp [(zv/Ba)(vw — v)] 


] 


Eq. 35 may be restated in the following theorem: 


Theorem 4. The velocity potential, in the mixed supersonic region on the top surface 
of a three-dimensional oscillating wing, may be computed by Eq. 35 or, in the z, y co- 
ordinates, 


U aoe 
&(z, y, + 0, 2) = —— exp (1) 
T 
e fs \ 2 .\ \2 2 2371/2 
x [| [Az, 7) + Ar, 1] cos L(y, 6 a)|( — £) = B (y — J") dn dt (35a) 
ES Q(a — ¢)? — By — n)]’” exp [(WU/6a)(z — &)] 


U , r A,(E, n) cos {(v/B’a)[(x — §)” — Bly — n)’|'?} dn dé 
- exp (iv?) | = acer aa | casas TF! ie [iso Vie 


(ec — & — By — 0)” exp [(U/6"a")(z — 8) ’ 


~ 
a 


Swe 
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where Sy, is the area bounded by 0 < u < w’ and v,(u) < v < »v,(u) and Sy, is the 
area bounded by u’ < u < uw and v,(u) < v < vy. 

The corresponding result for a point on the bottom surface of the wing may be 
obtained by interchanging A; and A, . It may be noted that for a wing of zero thickness, 
An(é, 7) + Arlt, 7) = 080 the integrals over Sw, vanish. The simple “equivalent area” 
theorem established by Evvard in Ref. 2 for steady flows is thus seen to be valid for 
oscillating flows. 

Discussion. (A) In Ref. 4, the boundary value problem for the determination of the 
velocity potential in the purely supersonic region of a wing in unsteady motion at super- 
sonic speed was treated by source-superposition method in a quite general manner. In 
fact, theorems 1 and 2 mentioned above are included in Garrick and Rubinow’s results. 
On specializing to considerations of an oscillating wing, the derivation of Eq. 16 becomes 
very simple. The derivation of the same equation in Ref. 4 is more complicated. 

(B) Since an arbitrary downwash function can usually be expanded as a Fourier 
series, a harmonically oscillatory motion may be considered as a basis for building up 
more general motion for a nonstationary wing. This paves the way to construct a proof 
for a theorem applicable to more general nonstationary wings. Evvard, in Ref. 1, treated 
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Fia. 5. Contour C in »-plane. 


the general mixed boundary value problem by the source superposition method. His 
results include the present theorems 3 and 4; however, his analysis was not carried 
through to the present point. 

(C) A particular type of motion which is of both theoretical and practical interest 
and which demonstrates simply that theorems 3 and 4 apparently do not apply in the 
simple equivalent area form to all nonsteady motions, is the so-called “unit step” 
motion, in which a wing at rest starts abruptly at a certain instant and then maintains 
a steady motion. For composition of the velocity potential for a wing with motion of 
this nature, the “unit step” source will be useful. The “unit step’ source can be derived 
from an oscillating source by a contour integration in the v plane, 
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a lf aw Bp f (27: - » U dv de 
$,(z, y, 2, t) = Oni | ay = Serr, | cos an) exp ‘= t{— Pa (a — &) 2 (36) 


where C is the contour shown in Fig. 5. By writing the cosine term in exponential form 


Eq. 36 can be shown to yield 


8B Ux , r; F Ux r; 22. 
g(x, y, 2, t) = ar, Ee ae a 2) -+ H(t — Fa Bal 36a) 


where H(y) is the “‘unit step” function having the properly that 


\! u> 0, 
H(u) = ' 
lo p<. 0. 


Now, draw a sphere of radius (at) enclosed in the circular cone from the ‘“‘unit step” 
source at (£, n, ¢), with the center of the sphere located at a distance (U2) from (£, 7, ¢) 


S> 








Fig. 6. Region of influence of the “‘unit step’’ source at (, n, ¢) at instant t. 


(Fig. 6). Then the region of influence of the source is divided into three regions (by 
Eq. 37), 


(1) In region I, the influence is equivalent to that of a steady source. 
(2) In region II, the influence is equivalent to that of a steady source of half strength. 
(3) In region III, no influence of the source will be felt. 


The region of dependence for a point (x, y, z) will consist of three similar regions. 
Consider the lift problem of a rectangular flat plate wing performing a “unit step” 
motion. Suppose that the velocity potential at a point S in the mixed supersonic region 
near the wing tip is to be computed at an instant ¢, , such that ai, < | y| (Fig. 7). Ac- 
cording to the above argument, the condition at S will depend on both regions A and 
B and the wing tip will have no influence. But in accordance with equivalent area form 
of Theorem 4 the domain of dependence at S would exclude the shaded region in Fig. 7 
in the computation of the velocity potential at S, at instant ¢, . This provides an example 
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of the fallacy of the equivalent area interpretation of Theorem 4 for arbitrary non- 
steady motions. 


Therefore, theorems 3 and 4 are not directly applicable to “unit step” wings. This 
fact is indicated (but not proved) by Eq. 36, because the operation of the contour inte- 
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Fic. 7. The wing tip region of a rectangular flat plate performing “unit step” motion at instant t . 


gration will carry the cosine function to infinity such that the argument of Eq. 25b 


breaks down in the proof of Theorem 3. 
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THE OSCILLATING RECTANGULAR AIRFOIL AT SUPERSONIC SPEEDS* 


BY 
JOHN W. MILES** 
University of California, Los Angeles 


Summary. The pressure distribution on a quarter infinite, zero thickness airfoil having 
a prescribed distribution of downwash (on the wing only), which exhibits a harmonic 
time dependence, is determined by a Fourier transform solution of the linearized, 
potential equation for supersonic flow. The solution is effected with the aid of the 
Wiener-Hopf technique and leads to a Green’s function, which may be expressed either 
as a finite, definite integral or as an expansion in powers of a dimensionless frequency 
parameter. It is shown that the results are applicable to the calculation of the forces 
and moments on rectangular airfoils of effective aspect ratio (A cot 0, where @ is the 
Mach angle) greater than unity. It appears that the force and moment coefficients of 
practical interest may be expressed in terms of known functions, including certain 
integrals which have been calculated for the two-dimensional, oscillating airfoil. The 
extension of the two-dimensional results to rectangular wings for which the prescribed 
downwash is constant along the span is particularly simple. The extension of the results 
for harmonic time dependence to the step function (Heaviside) case is indicated. 

1. Introduction. The linearized, two-dimensional problem of the oscillating airfoil at 
supersonic speeds has been studied by a number of analysts, using a variety of ap- 
proaches. This work has recently been collected and summarized in two reports prepared 
by Biot et. al.,'"* which give both the methods of analysis' and the numerical results.” 
(The reports also include the subsonic, compressible case.) These results are probably 
adequate for the strip theory analysis of wings with supersonic leading and trailing edges 
and aspect ratios sufficiently high to render tip effects small. (For an approximate treat- 
ment of the swept wing, reference may be made to a recent paper by the writer.*) Un- 
fortunately, those wings which meet the limitation of supersonic leading and trailing 
edges are generally characterized by small aspect ratio; moreover, the most serious 
supersonic flutter problems indicated by two dimensional analyses of such wings fre- 
quently occur in the near sonic regime, where tip effects are by no means negligible. 
Accordingly, it is of considerable practical importance to consider the three dimensional 
problem of the oscillating airfoil. 

The problem selected for study in the present paper is that of the rectangular wing 
tip, since it is the simplest three dimensional configuration (excepting those wings with 
no subsonic edges) of practical import. The results may be applied directly to the 
rectangular airfoil of aspect ratio sufficiently large to prevent the Mach waves from the 
leading edge wing tips from intersecting one another forward of the trailing edge and, 


*Received November 21, 1949. 
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2J. N. Karp, H. Weil, M. A. Biot, The oscillating airfoil in a compressible fluid, F-TR-1195-N D (1948). 
3J. W. Miles, Harmonic and transient motion of a swept wing in supersonic flow, J. Aero. Sci. 15, 


343-347 (1948) 
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indirectly, to the case where these Mach waves intersect on the wing but do not intersect 
the opposite side edges. It is also possible to solve the case of arbitrarily small aspect 
ratio, but the results are so complex as to be of dubious practical value. 

The method of solution to be used follows the Wiener-Hopf technique,* which has 
previously been applied to the rectangular wing in a steady flow.’ 

2. Statement of problem. A thin, rectangular airfoil is located in the vicinity of the 
plane z = 0, and its projection there occupies the first quadrant of the xy plane. A flow 
of supersonic velocity U is directed along the positive x axis, so that the leading edge 
of the airfoil is projected on the positive y axis and the port side edge on the positive x 
axis, as shown in Fig. 1. The boundary conditions are linearized in the usual manner, 
so that they may be applied at the projection of the airfoil on the plane z = 0, rather 


az 








Fic. 1. x, y, z axes and projection of airfoil on zy plane. 


than at the airfoil proper. The equations of flow will also be linearized,’ so that the 
problem may be subdivided into antisymmetric and symmetric cases with respect to 
the plane z = 0; only the former case is of interest here, since the latter situation does 
not give rise to lift. Accordingly, it is sufficient to consider the half space z > 0 and to 
apply the boundary conditions appropriate to the airfoil in the plane z = 0+. The 
problem to be solved is then the specification of the perturbation pressure over the first 
quadrant of the plane z = 0+ from a knowledge of the downwash there. 

‘N. Wiener and R. Paley, The Fourier transform in the complex domain, Amer. Ma. Soc. Colloq. 
Publ. 19, Ch. IV (1934); E. C. Titchmarsh, Theory of Fourier integrals, Oxford Press (1937) 339-349; 
E. Reissner, J. Ma. & Ph. 20, 219-223 (1941). 

5J. W. Miles, On the rectangular airfoil at supersonic speeds, No. Amer. Avia. Report AL 866 (1948). 

6A more complete discussion of the linearizing process and its various aspects is given by P. A. 
Lagerstrom, Linearized supersonic theory of conical wings, J.P.L. Progress Report 4-36. California Insti- 
tute of Technology (1947); NACA T.N. 1685 (1948). 
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The vector, perturbation velocity (q) due to the presence of the airfoil in the flow is 

specified as the gradient of a velocity potential Ub@, viz: 
q(x, y, z, 1) = UbV¢(a, y, z, 0) (2.1) 


(b is a characteristic length, to be chosen in any convenient manner.) The gage pressure 


follows from Newton’s law as: 


co ae 
p(x, y, z, t) = —pUb Di d(x, y, 2, b) (2.2) 
where D/Dtis the time differentiation operator in a fixed reference frame and, in linearized 
form, is given by: 
D 0 0 
—— on Cf a fi = 2. 
Dt Ox + él (2.3) 


The condition of continuity, after linearization, leads to the scalar Helmholtz equation 
in the fixed reference frame, viz.: 


1 D’ 
= 4 9 
o DE v1) (2.4) 


Vo 
where c is the sonic velocity for the ambient stream conditions. 

At this point, it is convenient to introduce the harmonic time dependence exp (iw), 
the Mach angle @, the frequency parameter «x, the dimensionless coordinates (2’, y’, 2’), 
and the modified potential, pressure, and downwash functions y(2’, y’, 2’), y(2’, y’), and 
a(x’, y’), in accordance with the relations 


6 = sin” (c/U) (2.5) 
xk = (wb/c) tan 0 (2.6) 
x = (b cot 6) x’ (2.7a) 
y = by’ (2.7b) 
z = be’ (2.7¢) 
d(x, y, z, t) = exp [i@t — xz’ ese 0—)]¥(2’, y’, 2’) (2.8) 
2 d(x, y, 0 +, t) = (U tan 6/b) exp [i@t — xx’ ese @)}y(2’, y’) (2.9a) 
0 ee 
y(2',y’) = (2, — ix sin ave’, y’, 0+) (2.9b) 
0 
a(x’, y’) = —35 We’, y’, OF) (2.10) 


Substitution of Eq. (2.8) in Eq. (2.4) yields the reduced equation 


Wy’ + Ware! _ Wrz! + ky (2.11) 
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while substituting Eq. (2.9a) in Eq. (2.2) yields 
p(x, y, O+, t) = —pU’ tan 6 exp [7(wt — xx’ ese 0)]}y(2’, y’) (2.12) 
The boundary value problem to be solved now may be posed as: find a solution to 
Eq. (2.11) which satisfies the boundary conditions 


¥.(2’, y’, O+) = —a,(z’, y’), y’>0 (2.13) 
y(x’, y’) = 0, y’ <0 (2.14) 
Y(2’, y’, 2’) = 0, Y tr (2.15) 


Eq. (2.13) states that the (modified) downwash a,(zx’, y’), is prescribed on the airfoil; 
Kq. (2.14) states that the pressure must vanish off the airfoil, since it is presumed asym- 
metric with respect to z, and only the airfoil is capable of supporting a discontinuity in 
pressure; and Eq. (2.15) states that the disturbance is propagated downstream and must 
vanish forward of the Mach waves originating at the leading edge of the airfoil. 

3. Fourier integral formulation. A general solution to Eq. (2.11) may be conveniently 
formulated in terms of Fourier integrals. The Fourier transformation of a function of 
the space coordinates (x, y) into its representative in the (u, v) spectrum will be denoted 
by a transition from lower to upper case letters in the functional notation, and the 
conjugate transform operators 7’ and 7' are defined by 
f(z, y) = T{F(u, »)} = a. | du | dvF(u, v) exp [i(ua + vy)] (3.1a) 


F(u, v) = T{f(z, y)} = = | dx | dy f(x, y) exp [—i(ua + vy)] (3.1b) 
In general, the parameters u and v may be allowed complex values, but the paths of 
integration in the uw and v planes must be suitably restricted in order to comply with 
both physical and mathematical requirements. Frequent reference will be made to 
Titchmarsh’ and Campbell and Foster’, simply by using the letters 7 and CF, followed 
by the appropriate equation number in the original source, although the notation used 
herein is not entirely consistent with either of these references.” 
In addition to the entire transforms of Eq. (3.1), it is expedient to introduce the 
notation 


Fi(u, v) = T{f(x, y)1(2, y)} (3.2a) 
F_(u, vy) = T{f(z, y)l(a, —y)} (3.2b) 


7E. C. Titchmarsh, Theory of Fourier integrals, Oxford Univ. Press (1937). 

8G. A. Campbell and R. M. Foster, Fourier integrals for practical applications, Bell Tel. Syst. Mono. 
B-584 (1942); also published by D. Van Nostrand and Co., New York, N. Y. (1948). 

In the case of transformation with respect to a single variable, we find it convenient, however, to 
use the notation of ref. 8, such that, e.g. f(z) = 7, {F(u)} = 1/20 f° du F(u) exp (iuz). 

Thus, the inversions of section 5 correspond to ref. 8. 
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I(x, y) = 1, z>0 and y>OdO 


(3.3) 


0, <= or y <0 


Eq. (3.3) defines the Heaviside step function in two variables, but if only one variable 
is indicated, e. g., I(x), the step is in that variable alone. F',(u, v) then represents the 
transform of a function which vanishes off the wing, while F_(u, v) represents the trans- 
form of a function which vanishes on the wing, recalling, cf. Eq. (2.15), that the solution 
vanishes identically forward of the wing. 

It is readily shown that an elementary solution to Eq. (2.11) is given by 


Yo(x, y, 2) = exp [t(ux + vy) — dz] (3.4) 


. [»? ro (u a °)]'” (3.5) 


the sign of the exponent Az being chosen to represent a disturbance which is bounded 
for large z. The primes have been dropped from the coordinates, but they are assumed 
to be the dimensionless (primed) coordinates of Eq. (2.7). By virtue of the linearization 
of the problem, these elementary solutions may be synthesized with the aid of the 
Fourier integral to form solutions capable of satisfying prescribed boundary conditions. 
In particular, the most general solution to Eq. (2.11) reducing to y(a, y), ef. Eqs. (2.9b) 


and (2.14), is given by 


er, PGs. s) exp (—W)} 
¥(z, Y; z) ee r i(u Pa x sin 6) (3.6) 


provided that the paths of integration in the u and vy planes are suitably chosen. 

In order to determine the appropriate paths in the complex transform planes, it is 
necessary to establish a domain of regularity for T,(u, v). If it is simply assumed that 
y(x, y) is bounded for large x (the behavior of y in y affects the behavior of T in v only), 
it follows from Eq. (3.1b) that I, (u, v) will be regular in Im (ux) < 0, so that the u inte- 
gration may be carried out along a path in the lower half of the complex yz plane. More 
specifically, I. (u, v) will be found to have a simple pole at 1» = +(« ese 6 + ze) if the 
pressure on the wing, ef. Eq. (2.9a) behaves as exp (— ex) for large x, where ¢ is a positive 
real constant, and will have a zero at ~« = « sin 6, as may be verified a posteriorz. (Due 
to the zero in T at » = « sin 6, V does not have a pole there.) In the present analysis, 
it suffices to take « = 0, insofar as Im (x) < 0. Accordingly, the integrand in Eq. (3.6) 
is regular in uw except for a simple pole on the real axis and the branch points in X, the 
latter being designated as +uo(v), where, cf. Eq. (3.5), 


h = (uo — w)'” (3.7) 


po(v) = (x° + "i led (3.8) 

The location of these branch points of course depends on », but if the path in the v plane 
is chosen such that 

| Im (v) | < « (3.9) 

they will always possess real parts, and their imaginary parts will always be less (in 


magnitude) than | Im (yv) |. It follows that, if the branch cuts from +: are both extended 
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to +7, the integrand of Eq. (3.6) will be regular everywhere in the half plane bounded 


by 

Im (u) <— | Im () | <— | Im (uw) | (3.10) 
and the phase of \ will be 2/2, 0, and —z/2 forIm (u) = O and Re (u) > wo , 
Re (u) | < wo and Re (u) < —ypo , respectively. An appropriate path of integration is 





then one satisfying the restriction (3.10), as illustrated in Fig. 2. 


nots Ru) 








= -& path for E¢. (3.6) | 





Centon pm: Re? J 
rrr x<2 Pa 
i ie. Ps 
ii: ae eee a 
Fic. 2. u« plane, showing branch points +o in A, the branch cuts therefrom, the path of integration for 
q. (3.6), and the completion of this path for x < 
Consider now the behavior of the solution (3.6) for x < z. The path of integration 


along the horizontal line chosen in accordance with Eq. (3.10) may be closed by a semi- 
circle in the lower half (x) plane. If u is denoted by Re’ along the latter path, the choice 
of phase for \ implic 


exp (7ux Nz ~ | exp [zul, z)] | ~ exp [—R(x — 2) sin 6] (3.11) 
Since @ is negative along the semi-circle, the integrand of Eq. (3.6) is exponentially 
bounded, as long as z, and the contribution of this path to the contour integral 
therefore vanishes uniformly as RF is allowed to become infinite (for —a + « < 0 < —e, 
e > 0; the portions —e < 6 < Oand —7z 6 << —a + «may be interpreted as belonging 


to the path associated with the transform path for y). Moreover, since the integrand is 


regular in the half plane bounded by Eq. (3.10), the entire contour integral vanishes by 
virtue of Cauchy’s theorem, whence y¥(x, y, z) vanishes identically for x < z, in satis- 
faction of Eq. (2.15). 

The asymptotic behavior of the solution for large, positive x (i.e., the Trefftz plane) 
may be obtained by closing the contour in Im (u) > 0 plus “keyhole” paths around the 
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branch points -+-w, , the only contributions to the integral coming from the latter paths 


and the pole at u = +« ese @."° 
Turning to the v plane, the branch points of \ may be designated as +») , where 


. = (v = pry”? (3.12) 


(a) = i — 2)! (3.13) 


Considering »)(u) in the » plane, the branch cuts are taken from +x to +, so that the 
phase of vy will be 0, —2/2, and —z for Im (uw) = 0 and Re (xu) > «, | Re (wu) | < «, 


[-K +K 7 Ri (Aa) 








M path (A) ‘é) (, ‘D) (E) 


Me par, 


(A (E) 
(B) ” 
©} 


Fic. 3. u plane, showing cuts from +x and paths of » and »(u). A, B, C, D, E, represent corresponding 





points on these paths. 


and Re (y) - x, respectively. As u follows the horizontal path designated by Eq. (3.10) 


| follow a lower path, as shown in Fig. 3, such that 
Im (v%) < Im (yu) (3.14) 


The cuts from the branch points + in the v plane must be chosen such that Re (A) > 0 
(in order that the solution (3.6) will be bounded for all positive z) for all (u, v); accord- 
ingly, the path of integration in the v plane must pass under one branch point and over 


WIf the path for Eq. (3.6) is allowed to approach the real axis (also taking Jm(v) = 0) and is indented 
under +x ese @ and po , the pole will give only half the contribution obtained for /m(u) < 0 and will 
also contribute an equal and opposite amount to the contour for x < z. It is necessary in this case to add 
an auxiliary solution, which is independent of x and is represented by ré(u — « ese 0)T.(x esc 8, v) 
* exp cot2 @)'/2z]. See, e.g., J. W. Miles, On linearized supersonic airfoil theory, No. Amer. 
Avia., Inc. Report AL-801 (1948) pp. 15, 16, where the care of stationary flow (x = 0) is discussed. 
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the other, and, in addition, the cuts must not interfer as Re (v)) changes sign, cf. Figs. 
3. and 4. It follows that the cuts must run from +» to —7@ and from —y, to +72, 
since Im (y)) < 0 for all u, as shown in Fig. 4. The result of this choice is that the ele- 
mentary solution Yo(2, y, z), ef. Eq. (3.4), behaves as exp [iu(a — z)] for large (absolute) 
values of u and not too large values of v, indicating the (Mach) wave front « = 2, repre- 
senting the locus of the waves originating at the leading edge of the wing.”’ 

Returning to the solution (3.6), differentiating with respect to z, and substituting 


in Eq. (2.13) yields 


a(x, y) = T{G(u, »)T.(u, v)}, y>0 3:19) 
G(u, vy) = —i(u — «sin 6)'[v — v(u)]'” (3.16) 
It may be remarked that Eq. (3.15) is valid for all y if a_(x, y) is added to a. (x, y). 


Eq. (3.15), for y > 0, and Eq. (2.14), for y < 0, together constitute a dual integral 


pIm tv) 
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Fic. 4. v plane, showing cuts from +» and the domains of regularity of the various transforms. 


equation (cf. 7, pp. 334-342) for the determination of y(z, y). Its reduction to transform 


equations may be effected by taking the inverse transforms of Eqs. (3.15), extended as 


valid for all y, and (2.14). The results are 
G(u, v)T.(u, v) = A.(u, v) + A_(y, v) (3.17a) 
I'_(u, vy) = 0 3.17b) 


liSee also Iq. (3.11) in this respect. 
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4. Solution of transform equation. The transform equation (3.17) may be solved with 
the aid of the Wiener-Hopf technique. While it is possible to treat an arbitrary distribu- 
tion a.(x, y), the difficulty of inverting the complicated transforms which arise in any 
general solution probably would render it impracticable. Accordingly, it will be assumed 
that 

a(x, y) = y” exp (—By)1(y)a, (zx) (4.1) 
which will be sufficiently general for most purposes. 

The transform of Eq. (4.1) with respect to y is given by CF 524.2, while the transform 
with respect to x will be denoted by A,(u),’” whence 


A(x, ») = (2m) "D(a + DB + 8)" AH) (4.2a) 
Antu) = Tylaal®)} = | df exp (—inb)au(® (4.2b) 


Accordingly, A.(u, v) is regular in v except for a pole of order n + 1 at v = 78. Moreover, 
it may be verified a posteriori that y(x, y) behaves as a,(x, y) for large y and, therefore, 
also has a pole of order n + 1 at v = 78, while the fact that y(z, y) vanishes for y = 0, 
ef. Eq. (2.14), implies that its transform must vanish at least as v™'~‘, where ¢ is a 
positive constant, for large v. Similarly, it may be verified that a_(x, y) has a singularity 
of at worst y ‘* for small (negative) y and vanishes uniformly for large y, so that its 
transform vanishes for both large and small y. It follows that sufficient conditions for 
regularity of the various transforms in the complex y» and v planes, and, therefore, suit- 


ably restricted domains in these planes, are 


A.(p, v): Im (u) <—8, Im (v) < B 
A_(u, v): Im (u) <—8, Im (v) > 0 
Pie. os Im (un) <—8, Im (v) < 8 


G(u, v): Im (u) <—8, | Im (v) | < | Im(») | > B 


The restriction Im (u) < —8 is in accordance with Eq. (3.10) and conveniently places 
the line Im (v) = —Im (vy) above Im (v) = 8, as shown in Eq. 4. It should be remarked 
that these conditions refer only to the individual transforms and do not necessarily 
imply regularity of the complete transform in Eq. (3.6), for which it is sufficient to 
require Eqs. (3.9) and (3.10). 

Since A. and I, are regular in the lower half v plane and A_ in the upper half, it is 
expedient to split G into two functions, which are regular, respectively, in these two 


domains; thus: 


G(u, v) = G.(u, »)/G_(u, ») (4.3a) 
G.(u, v) = —i(u — «sin 6)'(v + aa”, Im (v) <— Im (%) (4.3b) 
G_(u, v) = (v — ») 7”, Im (v) > Im (y%) (4.3¢) 


12See footnote 9 regarding notation. 
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It should be specifically noted that G, and G_ are not defined in accordance with Kq. 
(3.2). 

The domains of regularity of A, , A_ , T. , G, G, , and G_ in the v plane are shown in 
Fig. 4. Moreover, if Eq. 3.17a) is multiplied through by G_(u, v), it is evident that 


G_(u, v) A_(yu, v) is regular in Im (v) > 0 
G.(u, v)A.(p, v) is regular in Im (v) < B 
G_(u, v)A.(u, v) is regular in 0 < Im(v) < 8B 


Since A .(u, v) is regular throughout the entire v plane except for the pole of order (n + 1) 
at v = 78, it is possible to remove the singular part of G_(u, v)A.(u, v) to obtain a function 
which is regular in the upper half plane. Thus, introducing a Taylor expansion for 


G_(u, v) about v = 78, 


G_(u, v)A.(u, v) = | =. TC '(m + 1)G2""(p, 7B)(v — ia” Ja v) 


tL E u,v) — >> Cr '(m + 1)G*""(n, iB)(v — ia" | (u, v) 


it is found that the second term on the right is 0(1) at v = 7@, while the first term remains 
O[(v — i8)~"'], and it follows that the two terms are regular in the upper and lower 
half v planes, respectively. Multiplying Eq. (3.17a) through by G_(u, v) and subtracting 


from both sides the second term on the right of Eq. (4.4), the result may be written 


@,(v) = $_(v), 0< Im(v) < 8B (4.5) 


@.(v) = Gilu, vy) a(u, v) — | > r'(m + 1)G2""(u, 7B) — ia" |e v), 
“— (4.6) 


Im (v) < B 


@_(v) = G_(u, v)A_(p, v) + | Gl, v) 


_ i> T'(m + 1I)G2"(n, 7B) — ia)” | (u, v), Im (v) > 0 
where ®.(v) and ®_(v) are regular in the indicated domains, are identical in the common 
strip of regularity 0 < Im (v) < 8, and are, therefore, analytic continuations of the 
same function, say ®(v), in the lower and upper half planes, respectively. 

The problem is now reduced to the determination of the function (vy), and, since 
&(v) is regular for all v, it must be an integral function.’ Indeed, since it is regular at 
infinity as well as in the finite plane, it must be a constant, by virtue of Liouville’s 
theorem.'* To determine this constant, it suffices to examine ®(v) for large v. Recalling 


83. T. Whittaker and G. N. Watson, Modern analysis, Macmillan Co., New York (1947), p. 105. 
14T bid. 


~I 


1951] OSCILLATING RECTANGULAR AIRFOIL AT SUPERSONIC SPEEDS i 
that I, (u, v) must vanish at least as v™', that A.(u, v) is of order v™"~', and that G.(u, v) 
and G_(u, v) are O(v'’”) and O(v"'”’), respectively, it follows that #(v) is O(v"'”*) and 
therefore vanishes identically. 

Substituting the result 6 = 0, G, and G_ from Eq. (4.3), A, from Eq. (4.2), and 


evaluating the derivative of G_(u, v), viz. 
G""(u, v) = (—)"9 '?T(m + 1/2) — 7"? (4.8) 
the results for T, and A_ are given by Eqs. (4.6) and (4.7) as 


Po(u, v) = i(2r)'r 7 (u — «sin 0)(v + vm) (4B — wm)? (n + 1)A,(u) 


ee (4.9) 
1 . y— 6 » 2\—n=] 
- YI (m + 1)T(m + 1/2)( Bye + in 

aan 4 ice 1B 

A_(u, v) = (Qn)! '?'(v — )' 7 (G8 — v9) '?T(n + 1)A,(u) 
— ig\" (4.10) 

. > C'(m + 1)T(m + 1 2)( ~ 8) (8 + wv)" — A.(u, v) 
- oe! 


This completes the solution of the transform equation (3.17), and the remainder of the 
paper will be devoted primarily to the pressure distribution, as represented by IT’, (x, v). 
A similar treatment may be accorded A_(y, v) to obtain the downwash off the wing, 
but it is of less practical importance. 

In interpreting the result for I, (u, v), it is convenient to regroup the terms and 
write, cf. Eq. (2.9b), 


T.(u, v) = iu — «sin 6)¥(u, v) (4.11) 


W(u, v) = (Qn) er (iv + im) TE (n + 1A, (uy) 
: (4.12) 
- SS r'(m + YI(m + 1/2)(8 + in)" "(8 + iv)" 


where V(x, y) represents the potential only in the sense of Eq. (2.9b)."° 


In principle, the problem is now reduced to quadrature, i.e. the inversion of the 
transform W(x, v). Before considering this inversion, it is of interest to infer the behaviour 
of y(x, y) for small and large y directly from the behavior of I’, (u, v) for large and small 
vy, respectively. Thus, it is found that the pressure behaves as y'” for small y, since 
its transform behaves as v °”” for large v, while it behaves as y" exp (—y) for large y, 
since the transform behaves as (8 + iv)~"”* for small v, and the initial assumptions are 
verified. Similarly, it is found that a_(z, y) behaves as y~’”’ for small y and as 
((x + y)/—y)'” 1( 2, y) for large y, therefore vanishing outside the Mach line x = —y 
and behaving as (—y)~'” for large (negative) values of y within this Mach line, again 
verifying the initial assumptions. 

15Tnasmuch as the pressure distribution over the wing is independent of the discontinuity of y across 
the vortex sheet aft of the trailing edge (i.e., the “‘wake’’), it would have been possible to formulate the 
entire boundary value problem in terms of ¥. This would not be possible for a wing with a subsonic trailing 


edge 
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5. Pressure distribution. It appears that Eq. (4.12) cannot be inverted in finite terms 
of known functions. However, the result can be partially inverted so as to exhibit char- 
acteristic features, after which it will suffice to deal with certain integrals of the pressure. 

The » inversion of Eq. (4.12) may be effected by CF 524.2, CF 526, and the Faltung 


theorem. The end result is 


Wz, y) = 9 "'1y) > ( "\e(m + 1/2) [ dgala — 


Ss \m 
(5.1) 
| daly — 1)" kalE, 1) 
kn(é, 0) = T,{(iv)~"'” exp (—inoy)} (5.2) 


where (7,) is the binomial coefficient, and Eq. (5.2) utilizes the transform convention 
of footnote 9. Eq. (5.2) may be inverted with the aid of CF 571, CF 866, and the Faltung 


theorem, yielding 
: ee l 
k,.(z,y) =a T (2 = Nie, a y) | 


2 2 ii 


. 


dé|(x ve £)/2«)""*~* 4 
| (5.3) 


| 
a 
—~ 
ay 
—“ 
> 
w¥Y 


While Eq. (5.3) appears too complex for direct applications, it establishes the useful 
fact that k,,(x, y) vanishes for x < 0 or y > z. 

$v virtue of the fact that k,,(&, 7) vanishes for 7 > &, the upper limit y in Eq. (5.1) 
may be replaced by infinity in the region where y > x. Then, if (y — 7)" ” is expanded 
by the binomial theorem, and the 7 integration carried out prior to the uv inversion, the 
latter may be effected by CF 571 and the end result expressed in terms of powers of y 
and integrals of Bessel functions of integral and half integral order. This solution is, 
however, readily obtainable by more direct methods, since the region in question is not 
influenced by the region between the side edge and the Mach line y = —xz. 

For small values of y the » inversion of Eq. (4.12) may be brought about by ex- 
panding the transform in powers of (1/7), but the results would still involve Bessel 
functions of fractional (odd multiples of 1/4), therefore being only slightly less compli- 
cated than the result of Eq. (5.3). 

Still another approach is to expand the solution in powers of the frequency parameter 


x. Thus, if we determine the coefficients c,”’ in the expansion 


(ivy) '? exp (— ivy) 
~ . (5.4) 
* \—m=—1/2 . 2p (m) > -2p.a@ 
= (in)? exp (—iny) D> WK?” Do ck (in)? 
p=0 q=0 
the u inversion of Eq. (5.2) is given by CF 516 as 
. » | 
/ ” 2p (mm n-ll « _ 
k(x, y) = lly, zt — y) p> K by I'| 2p qtmt+ 5 , 
p=0 a=0 - er) 
5.5) 


2p—-q+m—1/2 


-y"(z — y) 
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If only the first term (p = 0) in this expansion is retained, the summation in Eq. (5.1) 
is binomial, and the result may be written 
z v : 
Vx, y) =21y,2—y) | ae l dno - 
0 


0 


(5.6) 
- a(x — & y + & — W)ME — my) + O(x’) 


in agreement with the result obtained in footnote reference 5. 
6. Wings of finite aspect ratio. Consider the rectangular planform depicted in Fig. 5, 











U 
Y. 
a ie es. 
67 ‘\ Fa eC: 
ok ‘\ fe 4 \e 
ts Xs 4 SP 
/ ar ° / \ 
No 7 
NN rd 
> 4 - & : 
\ 7 u N 2 
2 ! ~ Ps .t 
~~ & n 
x 
4 SN 
4 
4 Ss + 
/ x'=2tand 
Ps . . mn ey 
‘ 
a Po ™ \te 
4 vS 
1” Np 
$/ Ne 
yx’ _ 





Fic. 5. Rectangular planform and Mach waves. 


along with the Mach waves springing from the corners and the regions delineated by 
these waves. In the dimensionless coordinates of Eq. (2.7), the chord of the wing is 
chosen as 2 and the span as 2A. A’, the “effective’’ aspect ratio, is related to the true 


aspect ration (A) through 
A’ = A cot 6 (6.1) 


Due to the choice of the dimensionless coordinates, the Mach waves (or lines) make 
angles of 45° with the direction of the free stream flow. 

A point in region 0 is influenced only by points on the wing, since the fore Mach 
lines through such a point both intersect the leading edge, so that the downwash is 
known at every point in the zone of influence of such a point, and the pressure may then 
be determined by an integration of this downwash. The zone of influence of a point in 
region I, however, includes part of the region II, off the wing, where the downwash is 
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not directly prescribed, at least in the problem under consideration. The effect of this 
latter region is to supply a pressure deficiency in region I, which is just sufficient to 
cancel the excess pressure along the side edge y = 0. Accordingly, if c,,(z, y) denotes 
the pressure (jump) coefficient at (x, y) due to the downwash at points on the wing 
which lie in the fore Mach cone subtended from (x, y), and c,,(x, y) represents the pressure 
deficiency due to the points in region II within this Mach cone, the pressure in region 0 
is determined by ¢,,(z, y) alone, c,,(z, y) vanishing there, while the pressure in region | 
is determined by 

Cy(2, Y) = Cy(2, Y) — Cp,(2, Y) (6.2) 
Now, in the case of the quarter infinite wing, the only regions of interest are 0 and TI, 
and the required pressure coefficient indicated by Eq. (6.2) is given by Eq. (5.1), in 


accordance with the relation, cf. Eq. (2.12), 


pU~/2)'|p(x, y, O—) — pla, y, O+)] (6.3a) 


tan dy(x, y) exp [2(wt — xx ese 6)| 


That Eq. (6.3) does indeed exhibit a discontinuity (in its derivatives) across the Mach 
is indicated by the step function I(x — y), which occurs in the Green’s function 


The pressure in region I’ on a rectangular wing of finite span may be determined 


from the pressure in region I by invoking symmetry considerations. Thus, by virtue of 


the linearization of the problem, the downwash distribution over the wing may be 


broken down into symmetric and antisvmmetric distributions in accordance with the 


relations 


a re a r, 2A — y) 6.4a) 
a (v, y — a r, 2A — y) (6.4b) 
and the corresponding pressures determined by a*’(x, y) and a“ (x, y) will satisfy 


similar relations, whence it follows that 
Cc [s. Y) on tr. ¥ Fe 2. QA = y) (6.5) 


the top and bottom signs being associated with the symmetric and antisymmetric 


Accordingly, the pressure at any point on a rectangular wing, 


problems re specth ely 

for which the Mach lines from the leading edge corners do not intersect on the wing 

(i.e.. A’ > 2). or indeed for any wing which has the leading edge x = 0, the side edges 
- e > z P 


y = 0, 2A, and a trailing edge which lies entirely in the regions 0, I, and I’ and nowhere 


has a sweepback angle in excess of the Mach angle (45° in the dimensionless coordinates), 
g 
is directly determined by the solution for the quarter infinite wing. 


Consider, now, the less restricted case for which the Mach lines from the leading 


edge corners do not intersect the opposite side edges (i.e., A’ > 1). For this wing, the 
pressure coefficient in region III may be cast in the form 


c 2, y) = c¢,,(2, y) — ¢,,(z, y) Fe,(z,2A — y) (6.6) 


I 


This result follows directly from superposing the pressure deficiencies from the two 


edges. If the pressure coefficient given by Eq. (6.6) is now integrated over the wing 
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with a weighting factor which is restricted in its (2, y) dependence to have the same 
symmetry as a,(x, y), it is found that the result is identical with that which would have 
been obtained on the assumption A’ > 2. It follows that the results of this paper may 
be used for the calculation of the forces and moments on a rectangular wing subject 
to the restriction A’ > 1. 

7. Lift and moments. The integrals of primary interest for the oscillating wing are 
the lift, mid-chord pitching moment, and rolling moment. The appropriate, dimensionless 


coefficients are defined as follows: 


[ -1 e 
C, = E pUS exp (iwl) [| [—2p(x, y, 0O+)] dx dy (7.1) 


S 


eee rt on 2 7 
Cx E pUS-2b-exp (it | I| (b — x)[—2p(x, y, 0+)] dx dy (7.2) 

[y mas -1 pp 7 
c | 5 pPU'S-2bA- exp (iwt) I (bA — y)[—2p(x, y, O+)] dx dy (7.3) 

_ a ve 
here /) are the true coordinates, S is the wing area, 2b is the chord, and 2bA is 
the sp 

Whereas as more convenient in studying solutions to the potential equation to 


ie quantities defined in Eqs. (2.5)-(2.10), it is rather more natural in dealing 


ae il \ | 
th the forces on the rectangular wing to introduce the new, dimensionless quantities 
zo = 2/2b = 1/22’ cot 6 (7 4a) 
y } = y bh = y’ (7.4b) 
—_— » BD ae 
aX(x2, y®) exp (iwt) = —U™' Di 2 y, t) (7.5a) 
a(x2, y2) = exp (—”Axr?)a(2’, y’) (7.5b) 


a 4 ; " 
exp (wl £ pl ‘ [—2p(x, y, O+)] = 41 Di g(x, y, 0+, b) (7 .6a) 


y(x2, y2) = 4 tan 6 exp (—a”Axr7)7y(2’, y’) (7.6b) 
A = 2X see’ 6 (7.7) 
k = @b/U) (7.8) 


The unit of length selected is the half chord (6), in accordance with the generally ac- 
cepted, two-dimensional notation. The reduced frequency parameters k and ) are also 
defined in accordance with standard notation, namely that used by Biot."* \ should not 


loc cit. 1, 2. 
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be confused with the propagation constant defined by Eq. (3.5) and is henceforth used 
only in the sense of Eq. (7.7). a7 is the ratio of the (amplitude of the) downwash velocity 
to the free stream velocity, and y" is the dimensionless, pressure jump across the wing. 
In terms of this notation, the substantial time differentiation, cf. Eq. (2.3), is effected 
by the operator 
: a [ re] : ta 2s 
)= (l 2b) ( )+ Qk +) (7.9) 
Dt 820 
The dimensionless pressure due the downwash distribution 


: Yan(2”) 7.10) 


az, ¥ = (Y-) ay \a 


as 


is given by Eqs. (7 6b), (2.9b), and (5.4e) : 


a : 0 , ) MD, 
y(ar5, y?) = 4 tan Al(y-)r ( - + 2ik | ST )i (m + 1/2) 
} os \m 
7.11) 
| dé exp (—ar)e2(292 — 8) | dyn '?(y2 — n)""k,,(2E tan 8, n) 
The corresponding lift and moment coefficients, as given by Eqs. (7.1)-(7.3) and (7.6) 


may be written 


Ci, = | dx “¥2(2) (7.12) 


‘y(a) = A* | dy y'y2(a, y) (7.15) 


Substituting y2 from Eq. (7.11) in Eq. (7.15) and integrating by parts yields 


\ 


'y(a) = 41 (2 + 2ik > (")r + 1/2 
Yn\v) = 4Tan ie 2ik) DF (Mm + 4) 


(7.16) 


re 
i 


dé exp (—7Ab)ar(a — &)Kmn(2E tan 6) 


pA 


ko (2z) =a n—-mt+ a at | dy y "(1 —A 


y)""*"Km(a, y) (7.17) 


nA 


1 1 9\-! \—1 4n—m —1/2 
k(t) = 9 (n — m+ 2) '(n — m+ 1) A” | dy y 
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In order to evaluate the k’,,(x), the restriction A’ > 2 will be imposed (recalling, ef. 
section 6, that the end results for the force coefficients will be valid for A’ > 1), whence 
2x" tan @ < A. Accordingly, the upper limit A in the y integrations may be replaced by 
infinity, by virtue of the step function I(x — y) in k,,(x), ef. Eq. (5.3). Then, if k,,(x, y) 
is substituted in Eqs. (7.17) and (7.18), the y integrals may be evaluated in terms of 
Gamma functions and the u inversions effected by CF 571, whence 


n—m+1,2 


Rea) =a2'?Tn—-m+1) dD (-)'T'n —m—s4+ 2,3) (s + 1) 


a=0 


7.19 
1 1 x (m+s)/2 ( 
, fh m+sct+ » |r(s + )a--{Z) J cmaey/2( kX) 1(2) 
Substituting this result in Eq. (7.16) yields 
1/27 ) ; - 5 
yx) = 4 tan Ol(x)r'?T(n + (2 + 2ik) » r(m + 1); (m + 1) 
> (—)T'n—-m—s+24+ar(s+ pr | dm +s+ »| 
(7.20) 


1 \ - ar : r) (m+s)/2 
‘ r(s 5) A | dé exp (—arg(e sin 6) 


© es (XE sin Oa, (x — &), ¢=0, 1 


8. Downwash independent of ». In principle, any practical downwash distribution 
may be expanded in terms of the form y"a,(x), the only restriction being that the span- 
wise dependence be sufficiently continuous to allow a power series expansion to the 
desired accuracy. In practice, the evaluation of the terms in Eq. (7.20) and the integra- 
tions indicated by Eqs. (7.12)-(7.14) will be exceedingly cumbersome for large n. The 
simplest case is that of a downwash which exhibits no spanwise dependence, i.e. n = 0, 
which will be treated in this section. In this case only C, and Cy are of interest, so that 
only vy is required. 

Setting n = 0 in Eq. (7.20) and dispensing with the indices yields 


0 s i ' = 
yx) = 4 tan ana)(2 {- 2ik) | dé exp (—71Arb)a“(x — &£) 
OL “0 


(8.1) 
[Jo(AE sin 8) — (AA’ sin 6)~" sin (AE sin 8)] 


The first term in the square brackets corresponds to the two dimensional result,'’ and 
the second term may be interpreted as the correction for finite aspect ratio. Moreover, 
the two dimensional results (for lift and moment coefficients) may be expressed in terms 
of a set of integrals (required for n = 0, 1, 2, 3) of the form” 


1 
fn = frld, 0) = [ £" exp (—7A£).J (AE sin 8) dé (8.2) 
“J. W. Miles, The aerodynamic forces on an oscillating airfoil at supersonic speeds, J. Aero Sci. 14, 
351-359 (1947), Eqs. (9) and (12). 


Bloc. cit. 2, p. 2. 
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Accordingly, the two dimensional results are applicable to the rectangular wing (A cot 
6 > 1) if the f, are replaced by 
f(r, 0, A’) = fild, 0) — (AA’ sin 6)" f*F(A, 8) (8.3) 
fr(A, #—) = | £" exp (—7\E) sin (AE sin 8) dé 8.4) 
Integrating Eq. (8.4) twice by parts yields the convenient recursion formula 
\~ cos 6 f* = {[(A sin @) cos (Asin 8) + (7A — n) sin (Asin @)] exp (—2A) 
Ss.) 
— 6,(A sin 6)} 2nd fx, + n(n — 1) frr2 


As a simple example, and as an additional check on the results, the lift coefficient 


for a flat, rectangular wing at an angle of attack a, in a steady flow (sok = X = 0) will 
be calculated. Setting a a, and k = X\ = O in Eq. (8.1) yields 

y-(z) = (4a, tan 6)[1 — (7/A’)] (8.6) 
The corresponding lift coefficient is given by Eq. (7.12b) as 

Cy, = (4a tan @)|1 — 2A’) ] (8.7) 


in agreement with the well known result of Busemann.”” 

9. Moment due to roll. As an example of an antisymmetric problem, the moment due 
to roll will be calculated. If the angular velocity in roll about the midspan line is p, the 
dimensionless downwash distribution is given by 
(x, y) = 1(x, y)(pbA/U)(A™y — 1) (9.1) 


[The required pressure function for the determination of the rolling moment, cf. 
| 8 

+ 7 . . ir 1 . . . . . we oe 

Eq. (7.14b), is therefore — ‘yf + ‘y?/A; calculating this quantity from Eq. (7.23), 


substituting in Eq. (7.14b), and integrating the resulting terms in J,(Aé sin @) and cos 


(AE sin 0) by parts, the integrals may be evaluated in terms of the f, and f* of Eqs. 
(8.2) and (8.4) with the result 


Cs = aC’ O(pbA U 


= —2/3[(tan 6 + 2ik) fy — 2ikf,] + 2(AA’ sin 6)~'[(tan 6 + 22k) f# 


— Qik f*] — 2(\A’ sin 6)*[(tan 6 + 27k) fy 
+ 2k sec 6(2k ese 6 — 2i cos 6 — isec 6)f, — 4k sec 6 ese Of, (9.2) 


tan @ exp (—72A)Jo(A sin 6)] — (AA’ sin 6)" {2(tan 6 + 27k) f* 
+ 2ksec 6[4k esc 6(1 + cos 6) — 37 cos 6 — isec Ol ft 


— 4k’ sec 6 ese 6(1 + cos 6) f* — tan 6 exp (—2A) sin (Asin @)} 





9A. Busemann, Infinitesimale kegelige Uberschallstrémung, Jahrbuch der Luftfahrtforschung 7B, 


105-121 (1943). 
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10. Arbitrary time dependence. While the present paper has been focused on the case 
of harmonic time dependence, the results are, in principle, applicable to the calculation 
of the forces on the wing due to a downwash with an arbitrary time dependence, by 
virtue of the linearization and the well known Fourier theorem. 

An alternative approach to the treatment of arbitrary time dependence would be to 
use the response to a step function as the basic solution. In this case, the dimensionless 


downwash is presumed to be of the form 
a(x, y, t) = a(a, y)1[t — t(2, y)] (10.1) 


which is to say that the disturbance at (x, y) arises abruptly at ¢,(z, y). Suppose that 
the solution is placed in the form 
Ye, y,0 = | a f dnale — fy — 1, — WG, mleAG, 2) (10.2) 


t 


Ss 


where y~ is the dimensionless pressure jump. Suppose further that the harmonic time 


dependence problem is written 


a(x, y, t) = a“(z, y) exp (iwt) (10.3) 


y(a, y, ) = exp (iat) [ dt [ dn k(e — & y — 1, wa, 0) (10.4) 


Then, by virtue of the Fourier representation of the step function of Eq. (10.1), a and k 


are related by 
a(x, y, t) = (2r)~" lim | 


dw (€ + iw)"k(zx, y, w) exp (twt) (10.5) 


which is to say that the indicial admittance a is the Bromwich integral of the Green’s 
function k. Accordingly, the two approaches are complementary, and, having solved for 
k, as in the present paper, a follows from Eq. (10.5). The application to the gust loading 
of a rectangular wing has been given in a separate paper.” 

11. Numerical results. Numerical results for the lifts and moments on a rectangular 
airfoil due to plunging and pitching oscillations have been obtained and are available 
elsewhere.” 

12. Related papers (added in proof). The special case of section 7 has been treated 
independently by Stewart and Li** and by Stewartson.** The results presented in the 
latter paper are in agreement, with those presented herein, whereas those of the former 
are not. Moreover, the result (8.1) has been checked by the author, using still a fourth 


method.** Hence, it appears that the general method developed by Stewart and Li 


may be in error. 
2J. W. Miles, Transient loading of supersonic rectangular airfoils, J. Aero. Sci., 17, 647-652 (1950). 

21J. W. Miles and Irven Naiman, Aerodynamic derivatives for oscillating rectangular airfoils at super- 
sonic speeds, U.S.N.O.T.S. Tech. Memo RRB-32, Inyokern, Calif. (1949). 
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THE AERODYNAMICS OF SUPERSONIC BIPLANES* 


BY 
W. R. SEARS anv H. S. TAN 
Cornell University 


1. Introduction. The drag of supersonic wings increases rapidly with increasing thick- 
ness. This has led to some speculation about the potentialities of supersonic biplanes, 
which might afford structural strength and rigidity by virtue of their external structure 
and hence permit the use of thinner airfoil profiles than would be possible in a mono- 
plane. This brings to mind the possibilities, recognized for several years, of actually 
reducing the drag of wings by providing the proper wave interactions between the upper 
and lower wings of a biplane arrangement. That this can be done in the two-dimensional 
case, i.e., in a biplane of infinite span, was proved in 1935 by Busemann (Ref. 1), who 
showed that the drag (excluding viscous drag) can be made equal to zero for a biplane 
at zero lift. 

Clearly, it is of interest to study the aerodynamics of finite-span biplanes at super- 
sonic speeds, and especially to estimate the effects of the wing tips on the drag of a 
finite “‘Busemann biplane.’”’ In this paper we shall report briefly on an investigation 
(Ref. 8) of the aerodynamics of biplanes having rectangular wings of identical planform. 
To simplify the work, we shall use here the small-perturbation linear theory, in which 
all shock and expansion waves are replaced by Mach waves inclined at the free-stream 
Mach angle. Busemann, to be sure, did not make this approximation in his two-dimen- 
sional biplane studies; nevertheless, it should be permissible for the slender airfoils that 
are of greatest practical interest. 

In the linearized theory, the Busemann biplane arrangement becomes the one shown 
in Fig. 1, i.e., the top and bottom surfaces are flat, the leading-edge Mach wave of either 
wing intersects the other wing at mid-chord, and the airfoil slopes are related by the 
formulas, for x > c/2, 

Yi(x) = —Yi(x — c/2), H(z) = —YVi(x — c/2). 
The typical case is then simply that of two isosceles triangles pointing at each other. 

In this investigation, the Busemann relationship between gap, chord, and Mach 
angles shown in Fig. 1 will always be assumed, but it will not be necessary to specify 
the shape of the profile in deriving some general results. It will be shown that the velocity 
potential, including all interaction effects, can be calculated by means of integrations 
involving the wing surface slopes only. The general results will be applied to the numerical 
calculation of the wave drag, at zero lift, of the typical Busemann arrangement having 
triangular wing sections. 

2. Formulas for source distributions. The equation satisfied by the disturbance 
velocity potential ¢ in the linearized theory is 


B°b.2 — by —%2=0; 6 = M*—1. (1) 


where subscripts denote partial differentiation with respect to the rectangular Cartesian 


*Received March 1, 1950. 
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coordinates x, y, z. Here M denotes the free-stream Mach number, and the coordinate 
x is taken in the direction of the undisturbed stream. It has been assumed in deriving 
Eq. (1) that ¢, , ¢, , and ¢, are small compared to the stream speed, U. A consistent 
approximate formula for the pressure coefficient is 


C, = 2(p — po)/poU> = —29,/U, (2) 


where 7» , po are the pressure and density of the undisturbed stream. 


o\ = 


An elementary solution of Eq. (1) is the so-called supersonic source, $(2, y, 
[(x — £)? — B(y — n)” — Bz — §)*]), provided that the value zero is taken outside 
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Fic. 1. The Busemann biplane arrangement. 


of the Mach cone that originates at the point &, 7, ¢. For brevity, we shall adopt the 
following notation: 
u(z) = [(a@ — &” — By — »)’ — Be]. 

It is well known (Refs. 2, 3) that a continuous distribution q of these singularities over 
a surface parallel to the flow yields a solution satisfying Eq. (1) and the boundary condi- 
tion 0¢/dn = rq on the surface. Moreover, Evvard (Ref. 4) has shown how a distribution 
of these sources over a fictitious diaphragm at a wing tip can be used to account for the 
interaction of upper and lower surfaces of a monoplane wing. 

We shall adopt Evvard’s scheme here for the calculation of tip effects for both upper 
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and lower wings, placing a diaphragm at each wing tip and introducing the conditions 
that these diaphragms are stream surfaces of the flow. The potential at points on the 
top (7') and bottom (B) surfaces of the upper (wu) wing is given by 


+ 


d.7(2, y) = -| quru(0) dS, (3) 


é 


dual, y) = — | queu(0) dS — [ qiru(c) dS, (4) 


and there are analogous formulas for the lower (l) wing. The areas of integration S, on 
the wing under consideration, and S’ on the other wing, are shown in Fig. 2. 











-——— 


ove 


Fic. 2. Diagram showing areas of integration, S of wing considered, and S’ of other wing, including 


portions of diaphragms 


Now the integrations over portions of S and S’ can be simplified immediately by use 
of monoplane results. First of all, it is clear that, in all areas unaffected by biplane inter- 
action, the wing-surface boundary condition requires that g = Uo/z where a is the 
slope of the wing profile in the x direction. Moreover, Evvard has shown, that for mono- 
planes—and therefore for biplane regions unaffected by interwing interaction—the inte- 
gration over the diaphragm can be replaced by another integration over part of the 
wing. For any point forward of mid-chord, i.e., x < a, there can be no biplane interaction, 
hence it is convenient to write the relatively simple expressions for these points before 


going on to treat the interacting regions. 
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x < a: no biplane interaction: Here monoplane results are applicable. For both upper 
and lower wings, we have (cf. Ref. 4) 


o7(z, y) = am: | oru(0) dS — = (og — o7)u(0) AS, (5) 
WT sy “aT / S171. 
; — 

$;(2, y) = ms [ opu(0) dS — o | (or — op)u(O) dS. (6) 
W /sy a /S1r. 


x > a: We consider now a point on the upper wing, top surface. If the point lies forward 
of the Mach line from the tip mid-chord (outside of area N in Fig. 3), there is again 
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Fic. 3. Diagram defining notation used in calculations for the upper wing. 
£ £ PI g 
aa’b—N, PP’Ow—S;, wOP”’—Srz , 
uguo’'O—S,o, Ou*ui—S,’, Ou*v*—Sy,’. 


no biplane interference and Eqs. (5) and (6) apply. For a point in N, however, there 
exists an effect of the lower wing, transmitted through the interaction region of the tip 
diaphragm. We can write 


U tao. U ‘ 

d.7(2, y) = — o,7u(0) dS — — d,,u(0) dS. (7) 
“Sr WwW sy 

Here, and in subsequent formulas, we denote by A(é, 7) the slopes of the tip diaphragms 

of upper and lower wings. Thus, for any point on the top of the upper-wing diaphragm, 

gur is equal to Ud,,/7, and this value has been used in Eq. (7). In regions unaffected by 

biplane interaction (e.g. for § < a), A(&, 7) is the same as for a monoplane, and Evvard’s 
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results will be used for such regions. In interacting regions ) is still unknown, of course; 
its determination constitutes the main problem of this investigation. We shall postpone 
this to the next section, after writing an analogous formula for points on the bottom 
surface of the upper wing. 

All points of the bottom surface of the upper wing, for which x > a are affected by 
biplane interaction. Let Sj} and S}, denote the areas of the lower wing and its diaphragm 
that affect the point (x, y). The wing-surface boundary condition is 


U 9 ; l 
Qua(t, Y) + : | [ o,ru(c) dS + [ duc) as “= Oun(X, Y)- (8) 
“S'r “S'rr 


r 0c 


This is an explicit formula for q.,(z, y), involving only known quantities. It may 
be noted that in the region S{, , q:r has been put equal to UA, /r. Moreover, here X, is 
a monoplane value unaffected by biplane interference, and is therefore known from 


Evvard’s work. We now have 


d.3(2, y) == [ duee(0) dS 
“S7+Sy7r 


(9) ° 
- ZF am@as-2f rm as, 

w Js; © Jars 

where q,2 in S; is known from Eq. (8) and X, in Sj, is known from monoplane theory. 
Again the calculation of the diaphragm source distribution, g,, in S;; , is postponed to 
the next section. 

For the lower wing there are formulas exactly analogous to Eqs. (7), (8) and (9), 
which will not be written out here. 

3. Calculation of diaphragm distributions. The conditions that insure that the tip 
diaphragms will be stream surfaces are the conditions of equal slope and equal pressure 
on top and bottom. Since, as Evvard has pointed out (Rei. 5), the diaphragms of a 
rectangular wing tip are not vortex sheets, equal pressures imply equal values of ¢, the 
perturbation velocity potential. We have, then, in region §;; , 


Od 4) 
. On and or = ds. (10) 
Oz Oz 
The first of these equations leads to 

La 

— A(z, y) = Gur(x, y) 

T 

(11) 


U a 
-_ — qua(x, 9) -- “ lf. o,ru(c) dS + [ dn(c) as | 


r a 
The second Eq. (10) states that, in S;, , 


= | o,,7u(0) dS — c | d,,u(0) dS 
TW Jsy Wi sry 
(12) 


2 l 
. -[ g.ou(0) dS — = [. eo 5 rm(c) aS, 


T 




























72 W. R. SEARS AND H. 8. TAN [Vol. IX, No. 1 


where gz in S; and S;; is given by Eqs. (8) and (11), respectively. We have now an 
integral equation for the diaphragm slope X,, : for points x, y in S;, , 


. 
. 


2] ru(0)dS= | (oe. — ovr)u(0) dS 


’Srr ’ Sy 


. 


+ | orme)dS+] rm ds (13) 


S’7 


1 afr r 
—_ | u(O) = | | o ru(c) dS* + | Aw(c) as*| ds. 


I S’II 


There is an analogous equation for A, , which will not be written out. 

Eq. (13) is to be satisfied for all points x, y on the upper-wing diaphragm. For some 
areas, there is no biplane interaction, i.e., S; and S7,; vanish, so that the second and 
third integrals on the right side of Eq. (13) disappear. It is clear that for these points 
the third integral vanishes as well, since S,; and S,; do not contain any points é, 7 affected 
by interaction. Consequently, for non-interacting points x, y, Eq. (13) reduces to Evvard’s 
integral equation for the diaphragm slope of a monoplane (Ref. 4). 

4. Solution of the integral equation. Eq. (13) can be written in the form 


AwH(O) dS = F,(a, y) (14) 
“Sir 

for points x, y in S;; , where 2F,(x, y) denotes the entire right-hand side of Eq. (13), 

and involves only known functions. We now introduce the new coordinates u, v, measured 

along the two families of Mach lines on the wing in question: 


M M.,. 
u = 53 (€ + Bn), v= 9, — Bn), 
B I 


c= yet, we ye, 


a(é, 7) 28 (15) 
ei = 7 = _ 3) 
OU, V) M 
( g\ 2 2/ 2 2.2)-1/2 
uz) = i(z — &) — Bly — 9) — Bz} 
| ae ; reo 
= 98 {(u, — u)(v, — v) — M*2’/4}-'”. 
Our integral equation now takes the form 
ihe 1 “’* X(u, v) dv , 
| i 172 | uO 73 = Fu, , ) (16) 
Jo (4 — WwW Jy 1 — v0) 


for points u, , v; in S;;. 
The solution can now be found by means of the following process: 
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[ F(u; ,v;) du, _ [" dy [" _du___—sf” Atv) dv 
Jo (u' — u,)' ie (u’ — ae Jo (Wu — u)*’* Jd, (@- y)'” . 
ig du “* H(u, v,) du 
= 7 ee (17) 
Jo (U —UM) Jo (4 — Bw) 
say, 
[” Hw,» 4 [ ~ I ite, 0) @ 
= U, V,) du —xo vz s= ff 1, V,) du. 
Je J. (Ww — w)'?(u, — u)'” Jo = 
Differentiating this result with respect to u’, we have 
a " Flu, ,v,) du, ** X(u’, v) dv 
2, [ Mam a [| Mee (18) 
Ou Jo (U — &) : a (v, =a 


We now multiply both sides of Eq. (18) by (v’ — v,)~'”’, integrate with respect to 
v, , and exchange order of integration in a manner similar to that just employed. The 
result is 


[ mz, —_ ( : [ Fit 90) de) = [ A(u’, v) dv (19) 


172 \4 7 
vy — 2) Ou Jo (U — %&) 


which implies (dropping the primes) 


ior (2 [ F(u, -¥) di) dv, 


Ouso (U — Uy) 


Au, v) = 172+ (20) 


172 (v — »,) 

This solution can be used to calculate the slopes \, in regions of interaction. This 
completes Eq. (7) for ¢,7 , and, by use of Eq. (11), also completes Eq. (9) for ¢.2 . Eq. 
(20) constitutes a generalization of Evvard’s expression for the tip-diaphragm slope, to 
which, in fact, it immediately reduces when uw, v lie in a region free of biplane interaction. 

5. Calculation of the potential. Although the biplane problem is now completely 
solved in principle, the straightforward calculation of ¢, especially for regions of biplane 
interference, by substitution in Eqs. (7) and (9), is extremely tedious. Fortunately, as 
will now be shown, it is possible to eliminate entirely the integration involving \, in 


wr Ov J. 


these two formulas. 
In both Eqs. (7) and (9), the term involving \, is 


wk [Bt {” Mee 
MJ, (uy — wu)?’ J, (@ — 2) 


| r,2(0) dS = 73) (21) 


where now uw; , ?; lie in region S; . 
We return to Eq. (13), which holds for points in S;; , and write it in the form 


a lu . Mrz ,, 
| Gap Ow) = a7 ems), (Sn), (22) 
where 
G u, v;) = [ . du(u, ») sd — : [ Sun — Sue) 
Ju (v, — ») 2J-. © — 2) 
(23) 
1 3 __ dy Rs 
2U JL. (v, — v)'” ac® (u, 2) 
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and 
» USS ee . 
@(u,,%) = - Oo, 7TK\C) dS + A u(c) dS?. (2 t) 
Tv /S'] “S'rI 
Actually, ¢’ is the potential contributed at u, , v; by the lower wing. 
The solution of Eq. (22) can be written down immediately (Ref. 6); viz., 
Mo f* é'(u’,v,) du’ e 
Gu, »:) = =; | RUT (u Sv,). (25) 
2U ou (ui—- UU) 


Since Eq. (22) is correct only for points uw, , v, in S;,—i.e. for u, S v;—we must restrict 


u in Eq. (25) as indicated. 
Now for points outside of the interaction region, i.e., for u 


< M’°c’/Av, , the inter- 
v,) is zero. Thus G(u, v,) is also zero for u < 


action potential ¢’(w’, M*c"/4v, . 
We can now consider an integral involving G(u, 2); 1.e., 


s du 
i 3 Glu, 0;) == | aS 7 G(u, v,), 
(%, — &)~ 


~ M*c*/40,4 


[ du 


— 


\u, — Uw 


lie 


;) can be taken from Eq. (25): 


f rs) [ b’(u’, v,) du’ 


where kx S 
If x < v, also, G(u, 2 


M f* du 


| a ae P / fa Ba 
2 lu, — U LOU J 3920240, 


- ee 
(et: -= F (26) 


“ M*c*/404 


M f* ee fn a ae 
OT? ) \u ; V;) ae ; oe du 
2U J rc2/40 Uy — K k—U u, — U 
after some manipulation. Recalling the meaning of G(u, v,), (Eq. (23)), we can write 
Eq. (26) as 

v) dv 


. du aes r (1, 
ig 172 172 
df ae = & ‘“ vy, <= UY 


= M | o’ ic. v(* _ = ) | 4 _ | du’ (27) 
4 i K/ 2 u; — u 


2U «| M22 ~~ = 
] du {/° (on — o,7) dv l [ dv 0 
. i a ———755 — ¢'(u, v) ?. 
' 2 | (u, —u)”” (J. (vy, —v)” U Ji. WU y)'”* ac? 
Since the only restrictions on Eq. (27) are « S u, , and x S »,, it is exactly the result 


we need for Eq. (21), in which x = v7; S u,. 
We are now prepared to write complete expressions for the potential on top and 
bottom surfaces of the upper wing, by substitution in Eqs. (7) and (9). Let S;, be the 


portion of S; for which u S v, , as indicated in Fig. 3; then 


yf ! : Cw os 
d.r\zt, y) = -—— | a, 7h(Q) dS — _ | (tus = o,,7)4(0) dS 
T Js aT J Si 
Pe eer = 
—-=| u(0) — ¢’(u, v) dS (28) 
2m J S10+811 0c . 


av 1/2 
| re Vv, — uw’ 1 1 
+ - | ’(u’, v;) - — —— = == ————— 7 | du’, 
2 U, — V4; v.— U “, — U 


aT J M%c2/40, 
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WP . 
dus(t, y) = — c,2u(0) dS — : [ u(0) = ¢’(u, v) dS 
wT Jsy WJ s7+8r1 Oc 
Uf (un — our)u(0) dS + = | (0) 2 (u,v) dS 
+ Qr ! Sto r “on or 2r “Sro+Sit i ) 0c ? sa adie 
(29) 
Ge wa ti. = “| 1 1 ] . 
= 2r om ® (u ‘ o(% = “) Vv; = u’ oe _ = u’ du 


+ ¢'(zx, y). 

Formulas (28) and (29) permit the calculation of the potential, and consequently 
the pressure distribution, on the biplane. It is seen that, whereas we have succeeded in 
eliminating the integrals involving A, , for the upper wing, we are left with integrals 
involving \, , to be taken over certain interaction-free areas. In fact, if interplane inter- 
action of a higher order were encountered, such as an area of the lower wing influenced 
by interacting regions of the upper wing, it would always be possible to eliminate the 
\ integral expressing the last stage of tip interaction. 

5. Application of results: The wave drag of a finite Busemann biplane at zero lift. The 
general results obtained here have been applied to one typical practical case, to date. 
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Fic. 4. Spanwise distribution of wave drag near the tip of a rectangular supersonic biplane wing at 
M = ¥/2. The distance y is measured from the tip, and the chord is equal to 20. 


This is the case of the finite Busemann biplane, having isosceles-triangular profiles, at 
zero lift, at a Mach number of +/2. The computations have been carried out on a com- 
puting machine, and integrations have been done by planimeter. These are only pre- 
liminary results of an investigation that is still in progress. 
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These numerical results are shown in Fig. 4, where spanwise distribution of wave 
drag is plotted. The average drag coefficient of the biplane has been computed from 


Fig. 4 and found to be 


Total wave drag ' 
Cp = r= = 0.82387/A. (30) 
poU'sS 
where 6 denotes the leading-edge angle of the profiles, S the area of one wing, and A 
the aspect ratio of one wing. It is interesting to compare the wave drag coefficient of a 
rectangular wing of double-wedge profile, which is 


Cp = 4% (31) 


For the monoplane, 6 denotes the half-angle of the wedge. As would be expected, the 
ratio of biplane to monoplane wave drags diminishes with increasing aspect ratio. 

It is not difficult to show that the force coefficients calculated according to this 
theory for any Mach number 1, , can be extended to any other value of M by means 
of the following similarity rule 


Co(M) = Co(M,) (32) 
1—M 
The same correction would apply to the lift coefficient, C,(1/7). It is to be understood 
here that the coefficients Cp(M) and C,(M) do not apply to the same biplane as Cp>(M,) 
and C,(M,) but to a new configuration proportioned as in Fig. 1 at the Mach number M. 
In particular, the result of the present numerical work can be written 


Cp = 0.8235°/( AB’). (33) 
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STRUCTURAL ANALYSIS BY DISTRIBUTION OF DEFORMATION* 


BY 
C. V. KLOUCEK 


Prague 


The external loading of any structure produces, on the one hand, deformation and on 
the other hand, stresses which spread from the source of deformation in all directions 
offered by the given structure; the bending moment may be regarded as a function of 
stress. The deformation in any arbitrary form, e.g. joint rotation, elongation, etc., repre- 
sents in reality the single visible and measurable quantity. The expression of the external 
loading by means of fixed-end moments opened up the way for all the slope deflection 
solutions, as for instance, equations, balancing methods and the direct distribution of 
deformation. The balancing methods were based on the idea of alternately locking and 
unlocking the joints. The validity of both principles presented by Professors G. A. Maney 
and Hardy Cross follow from the fundamental law of super-position. 

In the course of the original derivation of the D. of D.-system (see the author’s earlier 
papers’ and books’), the direct moment distribution methods (Culman, Ritter, Suter 
etc.) were used as a pattern. Since these older computation systems are no longer cur- 
rently known, the author here presents the derivation of the two basic D. of D.-relations 
using newer methods which are first of all briefly characterised. 

Slope-deflection equations require first of all the establishment of all fixed-end 
moments ./''’* on the assumption that all the joints of the given structure are fixed. 
After the simultaneous releasing of all joints, the sum of all the adjacent fixed-end mo- 
ments acts as an impulse on the loaded joint. The deformation waves thus produced 
circulate through all joints until the complete settling-down of the loaded structure. 
The state of rest thus attained corresponds to the condition of equilibrium in every joint 
>> M = 0. If we substitute into this equation the basic relation between the bending 


moment 1/,_, and the joint rotations ¢, 
M,-. = £,-.(29, + Yn) + M ge (1) 


then the slope deflection equation for the location-fixed joint [s] and for n-adjacent joints 


[k] can be written in the form 
PsP + > (E,-n¢e) _ > is Mh yg (2) 


The notation introduced here is the same as that used in the author’s earlier publica- 
tions: 


£,-, = &- ,—<denotes the stiffness factor of the member [s] — [Kk] and has the value J/L. 
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p,—denotes the stiffness factor of the joint [s] and equals twice the sum of all 
member factors, p, = 27. e4 
6,—denotes the actual angular deformation of the joint [s]. 
y, = 2E£6,—denotes only a formal simplification for structures with constant modulus 
E. 

The balancing methods of angle changes by Maney, Goldberg, Cotten, Grinter, 
Morris, Kammuller, etc. constitute in principle the deformation counterpart to the 
successive moment distribution method of Hardy Cross. The balancing system e.g. for 
a structure with three elastic joints may be carried out in various ways e.g. successively 
[1] — [2], [2] — [3] or throughout the whole structure [1] — [2] — [3], and the like. If the 
chosen successive system has to be exact, it must express the influence of every joint on 
every other joint. By means of the sum of the partial angular deformation in any one 
joint [s], we obtain the resultant rotation ¢, , 

%2=¢t+¢e. +e + ete. (3) 
The magnitude of the individual increments ¢% , ¢;*, etc. is not a random quantity, but 
is governed by generally valid laws which can be expressed mathematically. 

The direct D. of D.-Method was originally worked out as the deformation counterpart 
to the direct moment distribution method by E. Suter.* The derivation could naturally 


3 








Fia. 1. 


be worked out from any of the other exactly valid systems, e.g. according to Clapeyron. 
Here the two basic relations, for the resultant deformation of the loaded joint and for the 
direct deformation distribution to the adjacent joint, will be obtained by an analysis of 


8Die Festpunkimethode, 1923. 
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the balancing principle. The individual deformation increments ¢; , ¢%", etc., which were 
originally obtained by successive distribution there and back, will be expressed as re- 
sultant sums, thus avoiding the balancing system. The derived relations will be checked 
by means of methods working with resultant deformations ¢, (e.g. by equations). 
Solution for two elastic joints. For the case where a single span of the given structure 
is loaded (Fig. 1), the initial supposition that every g = 0 leads to the determination of the 
fixed-end moment M{'*. The first stage of the balancing procedure, i.e. unlocking the 
loaded joint [1], is governed by the equation for one elastic joint p,g’ = M;"*; the first 
partial deformation of the loaded joint [1] is therefore, according to Fig. 1a, 
._ Mit 
os ea . (4) 
Pi 
The second stage (Fig. 1b) assumes the locking of the joint [1] in the deformed position 
for ¢; and the unlocking of joint [2]. The relation between the known, non-variable (fixed) 
value ¢; and the unknown ¢3 is given by the single equation for the second joint 
Pos + £29; = 0. The first partial deformation of the unloaded joint [2] is therefore 
z z E,-2 
ar = si (5) 
Po 
After a further alternation of locking we obtain by backward distribution the second 


deformation increment in the first joint 


z E,_9 
gi = —¢2*—> (6) 

Py 
Successive distribution (there and back) and addition of the corresponding partial values 


then gives the resultant deformations 


a=¢eate te t etc. . (7) 


anit += <— ~ ' (8) 
2 


From these well-known balancing systems we can easily derive the direct D. of D.- 
solution for the two elastic joints of Fig. 1. We express the second deformation increment 
>; of Eq. (6) in terms of the first deformation gj of Eq. (5), 


& 


rz zr ti» a / 
ea .. = wees (9) 
Pi Pe 


where the deformation constant a can arise only between two elastic joints, and its 
absolute value lies within the limits 0 and 0.25. The next deformation increment ¢;** can 


be expressed by means of Eq. (9) as 
g 
zz zz §$1-2 2 z 
1 = ¢) —~— = a-29) (10) 
and the resultant angle ¢g, at the loaded joint is obtained by substituting the partial 


deformations (9), (10), etc., into Eq. (7), 


i = ei(1 + a2 + Qj-2 + @i-s + etc.). (11) 








»~ 
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The deformation process between two elastic joints thus develops—in the same way as 
for moment distribution*—in the form of a simple geometrical series. On expressing the 


series (for a,;-.2 < 1) in the form 


(14+ Sai.) = — (12) 


i=1 —™ @j~2 
and substituting for ¢; from Eq. (4), we obtain the final expression for the loaded joint as 


in Fig. 1, 


ui 
a aed (13) 


The direct distribution of this resultant deformation ¢, to the unloaded joint [2] then has 
in this case the simplest form according to Eq. (8), 


£1 


= —p, 2. (14) 


Both the above expressions for g, and g, have been obtained by the addition of all the 
partial results from the completed balancing process [1] — [2]; the same expressions can 
be derived by elimination from the two slope-deflection equations. The sign of the 
individual angular deformations obtained by distribution, direct or successive, follows 
directly from the visual conception; joints which rotate in a clockwise direction are 
positive. 

Solution for three elastic joints. The determination of the exact effect of a third joint 
opens up the way for the solution of every straight-beam structure and hence we devote 
more attention to this problem. 

The above balancing system between two adjacent joints can be used successively 
also for three elastic joints as in Fig. 2. For the temporarily fixed joint [3], the balancing 


path [1]-[2] gives, according to Eqs. (12), (13) and (14), 


\= 








i=@ ; uv fix 
a =¢i 4) = — 1, (15) 
a(t T 2, ; :) pi(l —_ Q~2) - 
oe (16) 
P2 


For the fixed joint [1], the successive distribution [2] — [3] of the partial deformation ¢: 
gives the resultant Yo = Yo + 2 5 


¢2 = ei(1 +d ais), (17) 
i=1 


, from Eq. (16) gives 





which after simplification and introduction of ¢ 
E1—2 (18 
> =-o 2int 1 
ig ais p2(I = Qo) ) 


The resultant deformation of the third joint is obtained from the known relation 


a 
= -g, 2. (19) 


es. = Y2 
Pz 








4See e.g. Professor H. Cross, Trans. ASCE, 96, pp. 147, 150, 1932. 
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When using only the two successive distributions [1] — [2] and [2] — [3], the whole 
balancing process is not yet at an end; from Fig. 2 it is clear that for instance, the effect 


mfr 
wae (8 Ss 8 se 8 
OS Or @ Co-3 3) Or 





Mays 


ry = 



















2*% 


Fic. 2. 


of unlocking joint [3] on joint [1] has not been expressed, i.e. that all possible deformation 
paths according to Fig. 4 have not yet been pursued. Both described successive and direct 


© 





© 
ome Teen wee om, <. 


Fia. 3. 


solutions cannot exceed a maximum error of 2% for three elastic joints and are therefore 
well suited for the solution of practical problems; the incomplete balancing process 
naturally cannot be used for the determination of the exact effect of the third joint. 
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The creative idea of locking and unlocking the joints can be used in any arbitrary 
where the simplest form of the carry-over factor £/p is successively introduced 
*) _*) me 


sequence, 
The successive distribution, e.g. in the row 1—2—3—2—3— 


for the corresponding span. 


2—1 as in Fig. 3 can be expressed in terms of the constants a, 


Ss -sbe-ole-she-aBy-sbens | 
1 = ¢ ° ae = $10%1-2%2-3 « (20) 
P2 P3 P2 Ps P2P1 


In Fig. 4 there are indicated all the distribution paths for three elastic joints which 
0 + etc. The evaluation of the first 


give for example in the loaded joint g, = ¢i + ¢i 4 


Mite 
Uli) 1 dics g 52-3 Bi _y Sr= = 
“TO haan 


- §4- wv =f 








Fic, 4. 


deformation increment ¢;" is possible by only two paths 1—2—1 and 1—2—3—2—1, 


which, expressed in the sense of Eq. (20), gives 


¢1. = $1(ay~2 + Qt; 202-3). (21) 
Then the first rough approximation for g, = ¢; + ¢;° can be written in the form 
Mt R 
_i = ie jl + a,-2(1 + Q-3)}, (22) 
1 


which is in agreement with the exact expression (27) given below on substituting for both 
1. In a similar way we could obtain further partial increments in the joint [1] 
This procedure would be exact but ee 
2—1— 


indices 7 = 
and could derive the final expression for ¢, . 
ably time-consuming; the same may be said of the balancing system 1—2—3— 
etc. according to Table 1 

All the balancing systems for three elastic joints indicated here are only of schematic 
nature, since they take no account of the various intervals of time required by individual 
deformation waves. On the assumption that the balancing process is continued until the 
complete settling-down of the given structure (¢ = @), it is naturally immaterial 


1951 STRUCTURAL ANALYSIS BY DISTRIBUTION OF DEFORMATION 83 


what order we add together the individual deformation increments g = ¢* + ¢” + ete. 
While maintaining the basic condition that the effect of every joint on every other joint 





should be expressed, we can choose the balancing system shown in Fig. 5. 
| wv corry- carry- 

S ~ Over + & over “Y, 

7 2 3 
|= ¥, factor factor 








x F;- ie « 
eh Gt oe Re 
1 yr | £88 

2 


gy” ie Bs (4*Y™) 
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() 
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(4°) «fz 1. 9 
F2-3 


Ba) FG Te #5 








2. 3 
| ge a « 9-2 x ¢ Y;") 
F g 
gm aide ¥ x (gr gen) 2 
(y™) | ae — > _— ete. 
3. 























fix 


Making use of Maney’s principle (every ¢ = 0) we obtain the fixed-end moment 
and the first value gj. On using the simple carry-over factor £/p, the unlocking of joint [2] 


gives 
£ z £1~2 
i= ie (23) 
P2 
Completing the balancing of the value g; between joints [2] — [3] gives according to 


Iq. (17) the result g3(1 + >> a3_;). After substituting ¢g{ for ¢3 from the above equation 
(23), the fixing of joint [2] and unfixing of joint [1] gives the first increment in the loaded 


\ ae 


joint 


zz my —2 i Fo _ 
Yi = #1 ao (1 + eae ~y 


(24) 
= ¢i%-2(1 + ) Q2~3). 
In the same way the second and third cycles give 
¢i = ¢1 %-2(1 + p 3-3) 
(25) 


1 [a -2(1 + > a3-s)]’ 


gr = gila:(1 + >> a3-s)]’. (26) 
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The exact value of the resultant deformation follows from the sum ¢, = ¢; + ¢; + ete. 


Pi 


according to Fig. 5, 


The balancing process for the case of three elastic joints thus obeys the law of a geometri- 


My" 
a_i SG f-2 iy £2-3 % —f (NV) 


A LA Qy-2 A Qo.3 Ls jy? 











ge 





B+ Sai) 7 


esc. 


Fore. Bop NPEas 4) 








cal series of higher order i.e. a series of a series. After expressing the inner series in finite 
form, the outer series can be expressed in the form 


— a 7 ‘ 1 
od — : (28 
T =| | . = Q\-2 (l= Qo—3) ) 


| ei 
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The exact value for the loaded joint is thus given by the simple expression 





A 


¥1 


M;"* 
= ! sa (29) 
pill = @-2/(l — 2-3) | 
By addition of the partial values according to Fig. 5, we obtain the exact result for the 


distributed deformations 





£2 ; 
= -—y -, (30 
" p(1 _ Q2—2) ) 
fo-3 ‘ 
sx = G2 ’ (31) 


Ps 


which differ from the previous expressions (18) and (19) only by the exact value ¢, . 
The results (29) to (31) can easily be checked by any other method derived from the 


final stage of the balancing process, i.e. from the condition of equilibrium > M = 0. 
Using the elimination method’ of Professor J. B. Wilbur, we have for M3-1vy = —M3-2, 


according to Fig. 5, 
£3-1v(2¢s) + £,_3(2¢; + 2) = 0 
32) 
, ee (32) 


ye = 
, 2(E5_3 + £3_1v) ? Ps 


The condition of equilibrium >> M, = 0 for the joint [2], expressed in terms of deforma- 


tions 





£,-2(2¢, + $1) + f-a( 20, — $2 2) = 0 
3 


. . ‘a 2 
gives, after rearrangement and for £-3/pop3 = Qe-3 ; 


é 
—i (33) 


ee eS : 
pl — Qp-3) 


The condition of equilibrium 1/,_; + M,-. = M‘™ expressed in terms of deformations 


£,-1(2¢,) + ff Do — ¢ — fs — | = mM‘ 
p2(1 a 2-3) 


gives, after introducing the constant a,_, , the known expression (29), 





mi 
ig ie . (34) 

pill — a~2/(1 — a-s)] 

If for the structure of Fig. 5, we introduce an elastic support at the fourth point 
(denoted [4] instead of [IV]), the deformation of the loaded joint [1] proceeds according 
to a geometrical series of the third order. Successive introduction of the finite limits for 
the individual series results in a continued (chain) fraction, which we denote by the chain 
value a’, 
ut M?'** 

a; wes ! . . (35) 
pi(l — Q~2) pi {1 — @)~-2/ {1 — Qe-3/ (1 ae a3—4)]} 





’Transactions A.S8.C.E., 102, p. 346 (1937). 
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In the same way the first deformation distribution to the joint [2] will be 


tr 


£1-2 


ia p2(1 = ao-3) i, po[1 —™ 2-3 (1 = 3-4) } y 


1-2 





$2 


(36) 


For the case where the central joint [2] is loaded, we can evaluate the exact value ¢, 
by the sum of two balancing systems A and B as shown in Fig. 6. The simple distribution 
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Fic. 6. 


[2] — [1] — [2] of the value g3 gives 93° = ¢g}a,-. and after completing the balancing 
[2] — [3] of the value ¢3*, we have the first increase 


g2 (1 + bm a3) = g3a,-2(1 + 7, Gal, 
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and similarly 
On interchanging the subscripts, we obtain an analogical expression for the balancing 


g:"(1 + >> ai_s) = gilair-2(1 + >> a-s)]’. 


system B. The resultant deformation g, is obtained by the summation of two geometrical 


wot + > ai.) | (37) 


° we obtain as the 


seri S ol second order. 
(14+ $a.) [+E 
1=] t=! 


¢ = ef) | IE 


write all the series in finite form and rearrange the terms, 





When we \ 
deformation of the central joint [2] the simple expression 
M* 
= remmeenetarane ty (38) 
p(l — ay-: Q@2~3) 
This relation can be easily checked by means of slope-deflection equations. If we 
substitute for the two unknowns from the first and third equations 
E12 
G1 = 92° 
Dr 
(39) 
- S2-3 
ee 
Dea 
. . - - fix . . 
into the second equation prog, + £-2¢, + &-3¢3; = M,.'*, the latter may be written in the 
Mi 
sur sar pene. te (40) 


ale; a ae ge 
' pl — &:-2/ pipe — &2-3 P2 Ps) 


form 
For more complex open structures such as that shown in Fig. 7, we obtain by elimina- 

















~ 
Fia. 7. 


gtix 
M, (41) 


tion, in the manner of Eqs. (39) and (40), the expression for the loaded central joint [2], 
a) = 
7 poll — oie)’ 





®See Distribution of Deformation, pp. 26, 31, 32. 
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which is identical with the approximate formula derived by Mr. C. A. Willson, M.ASCE 
for a frame network.’ 
For the loading of the outer joint [3] of Fig. 7, we have 





M;" O39 
ee where a. = 3-2 ' (42) 
Ll = Gigg = Gna ~~ Gens 


p3(1 — a2) 


and for the fixed joint [3], i.e. for a2; = 0, the distribution of the known (non-variable) 


i 


deformation ¢g,; to the joint [2] gives 


r 53-2 (43) 


De = — S)- , . 
si ™ pol — ae-1 — 2-4 — 2-5) 





If we summarise now the results of the foregoing analysis, we see that every straight- 
beam structure can be solved by means of two generally valid relations for the loaded 
joint [s], 


>» M? x 





a (44) 
Au = b a,) 
and for the direct (resultant) deformation distribution to the adjacent joint [k]. 
_— Seat (45) 





a ea ae 
where n denotes the number of elastic members meeting in the joint under consideration, 

In concluding this paper, the author would like to add the following few remarks: 

This article has been devoted only to the theoretical analysis of the basic relations. 
The solutions for simultaneous loading, closed structures, joint displacement, etc. are 
given in earlier publications.” 

The author’s attention has been drawn several times to the fact that the solution of 
linear slope-deflection equations by means of continued fractions is of certain interest 
from a purely mathematical point of view; the author has not gone into this problem in 
detail but doubts if the mathematical applicability of the D of D.-system will be so 
general as, for instance, the Relaxation Method by Professor R. Southwell. 

The limited English edition of the work ‘‘Distribution of Deformation” has been 
printed only for specialists abroad; several copies of the book are in various American 
libraries (e.g. Brown University). The author intends to use the valuable information ob- 
tained through contact with specialists in various countries either for the working out of 
a comparative study or as an addition for the standard edition of the D. of D.-book 
where also the appropriate literature will be listed. 


7Transactions A.S.C.E., 102, pp. 352-354 (1937). 
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—NOTES— 


ROTATION OF AN INFINITE PLANE LAMINA: BOUNDARY LAYER GROWTH: 
MOTION STARTED IMPULSIVELY FROM REST* 


By SWAMI DAYAL NIGAM (Agra College, Agra, India) 


1. Introduction. T. v. KAirmén‘ has solved the problem of rotation of an infinite plane 
lamina in a viscous fluid. He assumes that the motion is steady and the lamina rotates 
with a constant angular velocity 2 about the axis r = 0. He has found exact solutions of 
the equations of motion which satisfy all the boundary conditions of the problem. The 
axial velocity does not vanish at infinity, but tends to a finite negative limit, which signi- 
fies a steady axial flow towards the rotating lamina. v. Karman interprets that it is 
necessary to preserve continuity, since the rotating lamina acts like a centrifugal fan, 
the fluid moving radially outwards, especially near the lamina. 

In the present note I have discussed the growth of motion in the earlier stages of its 
development caused by an infinite plane lamina which at t = 0 is suddenly made to rotate 
with a constant angular spin 2 about the axis r = 0. There grows a boundary layer of 
thickness proportional to the square root of time, adjacent to the rotating lamina which 
initially has a zero thickness. 

We start with the equations of motion in cylindrical coordinates and substitute in 
them expressions for u, v, w and p somewhat similar to those used by v. Karman. Then 
applying the approximations of the boundary layer theory, we integrate them analytically, 
satisfying all the boundary conditions required by the problem. The solutions have a 
serious limitation in that they give initial motion only. They give no information regard- 
ing the time after which the steady state is reached. 

2. Equations of motion. The equations of motion in cylindrical coordinates with 
terms of azimuthal variation omitted are 





au, lou u, &u ldp du du du wv 
Js tee He tes er Se Se te t+ 
| or” + r or r a] p or ot , or + Oz 7 
Ey 1 ov v av Ov Ov ov , w 
2 -—— = —— = ad eo <-- — 
” ar 7; or 7 a4 ot sis or + Oz r (1) 
jaw, 1 ow rs _ lop _ dw, dw, | dw 
| Or r or Oz p Oz ot or Oz 





where u, v, w and p are the radial, azimuthal, and axial components of velocity and press- 
ure respectively. The equation of continuity is 


re) 0 
> (ru) + a (rw) = 0 (2) 


*Received Oct. 14, 1949. 
1Goldstein, Modern developments in fluid dynamics, vol. I, The Clarendon Press, Oxford, 1938, p. 111. 
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We now let 


u = 2 rtf(n) 


v = Qrg(n) 
(3) 


w= —47' 7? h(n) 


p = 2wpQ*tp(n) 


where 
- 1/2 
n = 2/2(vt) (4) 
During early stages of motion when ¢ is small (or in boundary layer theory terminol- 
ogy: when the thickness of the boundary layer is small), we may neglect the terms in the 
equations of motion containing higher orders of t. Therefore by omitting terms of order ¢? 
in the equations of motion and continuity, we get to a first order of approximation, the 


following equations” 


fl 4 nf! — 4f = —4g? = —p” (5) 
g’”’ + 2ng’ = 0 (6) 
h’’ + 2nh’ — 6h = —p (7) 
f=h’ (8) 


3. Solutions of the equations. From (6) we get 


g = [1 —erfn] =erfen (9) 
With this solution for g the general solution of Eq. (5) may be expressed® 
f= A(1+27)+Bi(i + 2n’er fe : oo Qa ne “] + Wa *e" — n er fe n) (10) 


The boundary condition that f = 0 at 7 = 0 and 7 = & gives 

A = 6, B = 2/s 
The function h is obtained by a quadrature of f, and the function p by the double quadra- 
ture of 4g’. The final analytic expressions for f, h and p, are 


= - , ? 7 2 —* ¢ -—1/2 —y* o 

f = [(1 + 2n )erfe n — 2x ne) — Wr '7e" — nerfe ) (11) 
= 2 ; ; —1/2 | 2\ —7? 2n ne a. 

h = = [8n + 2n)erfe n — Qe (1 + ne" | — 7° <= n er fe ) 


2 oe : 2 2 2 
— ae Verfe n + xia erfe 2'"'n + aaa (= — 27 + 1). 


2/..\1/72 
3(7) T 


2Goldstein, loc. cit. p. 183. A similar approximation has been made there. 
The author is indebted to the referee of this paper for the particular solution to Eq. (5). 
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1/2 
p = (1 + 2n’\(erfe n) — = e" erfe n — a ow + 4(2) n er fe 24-1 + 2) 
T 7 7 T 
(13) 
2 /9 
oh (2 — 2? 4 1) + const. 
T vs 


Note that there is an anomaly in the behaviour of the pressure which makes it approach 
infinity as 7 —@ and precludes the specification of the constant in the last quadrature 
by a reasonable physical boundary condition. The anomaly is due to the acceleration of 
the infinite mass of fluid. 

The functions g, f, and h are given in Fig. 1. 
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Fig. 1. The flow functions. 


4. Stream function. A stream function may be defined by the equations 


ns 
aha r dz’ 
(14) 
1 oy 
ars 


whence the stream function may be expressed 
1/2 g2 43/2_2 ‘ 
Y= — 2 ONE’ h(n). 


This gives stream surfaces which are surfaces of revolution. 
In conclusion I express a deep sense of gratitude to Prof. M. Ray, D.Sc. for his kind 
help in preparation of this note. 
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RELATION BETWEEN BERGMAN’S AND CHAPLYGIN’S METHODS 
OF SOLVING THE HODOGRAPH EQUATION* 


By T. M. CHERRY (University of Melbourne) 


When a perfect gas is in steady irrotational isentropic motion in two dimensions, the 
stream function y satisfies a linear differential equation in which the independent vari- 
ables are components of velocity. For this ‘hodograph equation’, general forms of solution 
have been given by Chaplygin’ and Bergman’. The purpose of this note is to show how 
Bergman’s form of solution can be converted into Chaplygin’s. Hereby we obtain the 
specification of the same solution by means of two quite different series, and are in the 
position to check the extensive computations which are required (in general) to evaluate 
either of the series. 

The results of §1 are due to Chaplygin’, Lighthill* and Cherry‘; for proofs of the key- 
formulae (4), (6), (12) reference may be made to [3] or [4]. For Bergman’s form of solution 
the most convenient reference is v. Mises and Schiffer.’ The different authors use different 
notations, and the present paper uses a blend of them. 

1. Let the rectangular velocity-components be 7r’” cos 6, r'” sin 0, with the unit of 
speed so chosen that the limiting speed, at which the pressure vanishes, corresponds to 
7 = 1. Then the hodograph equation is 

42 2 2 
ay ov) , Ar oy fr ~ _ ar) ZY 0, (1) 


sa — (7 547% ~~ ihe y¥— 1/ 06 





where y is the adiabatic index of the gas. This equation is soluble by separation of the 
variables, leading to Chaplygin’s form of solution 


v= DewvA(re™, (c, constant) (2) 


’ 
where v can take any real value except —2, —3, --- , and 


v,(7) = 7’ F(a, ’ b, 3V + 7), 





1 v(v + 1) 
a,+b,=v—- ) SS eee 
re ary 2y — 1)’ 
F denoting the hypergeometric series. 
For 7 fixed, ¥,(7) is a meromorphic function of »; its poles are at » = —2, —3, ; 
and its residue at vy = —m is —h,,y,,(7), where 
I'(a,,) (1 + m — bn) 
(3) 





hn = (a, — m)T(1 — b,,)T(m)TU + m) 
*Received May 15, 1950. 
1§, A. Chaplygin, Sci. Ann. Univ. Moscow, Phys-Math. Div. Pub. No. 21 (1904). 
2S. Bergman, N.A.C.A. Tech. Note No. 972 (1945). 
8M. J. Lighthill, Proc. Roy. Soc. London, (A) 191, 342 and 352 (1947). 
4T,. M. Cherry, Proc. Roy. Soc. London, (A) 192, 45 (1947). 
5R. v. Mises and M. Schiffer, Advances in Applied Mechanics, vol. I, Academic Press, Inc., New York, 


1948, p. 249. 
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For large values of | v |, ¥,(r) changes character at the point tr = 7, = (y — 1)/(v + I) 
at which the coefficient of d°y/d@ in (1) vanishes; it is (for v real) monotonic for 
0 <7 <-+,, oscillating for 7, < + < 1. In the first case (with which alone we shall be 


concerned) there is an asymptotic expansion 
¥(7) ~ V(r) de af a p> parr}, (4) 


valid for all complex »v except the negative integers; here 


ee 1/2 — 1/2 
) — arc tanh (L =f) : 
T l-—r 





\ = 7,” arc tanh (z: _ 
(1 —_ aerennces a" (5) 
WI ; is ——-, Qa = 7,'” — 1, 
(1 — 7/7.) (1 + a)'** 


and the p,(r) are determinate functions vanishing at t = 0. Hence follows the partial 





V(r) = 


fraction expansion 


a vy wr = hate #a(0)| 
v.(r) = oe {%() ate” ee (6) 
valid for 0 < 7+ < 7, and all v. If here we formally expand the last factor in powers of »™* 


we must obtain (4), and therefore 
(—1)” Pn (7) V(r) = 3 m”~ Tm 5"e Yn (r), (n - 1, 2, r -*). (7) 


Substitute (6) in (2) and interchange the order of the double summation. We obtain, 
for 0 < T < Ts ? 


»_v(A+ 60) = — c,é- 
= V(r) > c,d’e ~ a hin d”e Wn 7) > mae 


C,  ) athdaaiieaad = c,8” 


= Do had eo Yale) rare 


m=2 





Hence, putting 


v(r, 6) = vhs sey .s) , (8) 
f= A+ 18, (9) 
Gols) = Di c,d'e”, (10) 


we obtain 


Vth 


= f 
V(r) > Cc, 5’e Oe oe »» hnoe™ Wn(7) [ y amnrr-= dt 


* t 
V(r)b0(8) — Do had" Yea) [ en" '-Pg (8) dt. (11) 


For the justification of these manipulations it is sufficient—apart from the over-riding 
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condition 0 < + < 7,—that (i) for all values of »y comprised in (2) and all positive integers 
svt | 


m, |v -+ m| has a positive lower bound, and (ii) the series >> | c,6’e” | converges in some 
strip ¢; < Re ¢ < 0; the proof rests essentially upon the estimations 


v,(r) = V(r)e’e"{1 + OW" )}, 2th, = 5°"{1 + O(m™)}, (12) 


of which the former is the first approximation derived from (4). 
We note that y,(7, @), as defined in (8), is a solution of (1) in Chaplygin’s form. 








2. To convert (11) into Bergman’s form of solution we expand the factor e”“~” 
and rearrange the resulting double sum. After an appeal to (7) this gives 
is r 4 
Y + vi = V(r)b0(5) — Dy mde” Yn(7) [ go(t) dt >) m™(t — 9)""/(m — 1! 
2 0 1 
t @ cs ; 
= V(r)¢(f) — [ > (t — §)""do(t) dt/(n — 1)! DS m""h,8"e" y,.(7) 
/9 n=1 m=2 
. , Ae — t"p,(1) V2) 
= Virbo(t) + | dS PRAT o(t) dt (13) 
#0 n=1 (n _ 1)! 
The transformation is valid provided the series 
z. h,o my (rent 
converges absolutely, and by (12) this is secured if | ¢ | 2 is negative; hence from (9), 
it is sufficient that \ be negative (as it is for 0 < 7 < 7,) and that 
—~"liaj<ce@ecs” iat. (14) 


On the left of (13) Y + y, is a solution of Chaplygin’s form, and on the right we have 
this expressed in Bergman’s form* in terms of an arbitrary analytic function ¢)(¢). The 
identification not merely of form but of content will be complete provided Bergman’s 
G, and the present p, are related by 
G, = (—2)” p,. (15) 


> 


Now if, as in [5], we examine the conditions that the form on the right of (13) be a solu- 
tion, with ¢,(¢) remaining arbitrary, we find that the derivative of p, must be determined 
entirely by p,-; , 80 that p, is determined apart from an additive constant. This constant 
is, in'the preceding work, determined by the condition p, = 0 for r = 0, while in [5] the 
condition is taken to be G, = 0 for \ = — ~; and these conditions agree since to r = 0 
corresponds by (5) \ =— ©. Hence (15) expresses merely the same function in two 
different notations. 

In conclusion, it may be remarked that the conditions assumed in proving (13) are, 
in one respect, more restrictive than those which validate Bergman’s form of solution on 
the right; for our conditions imply that ¢(¢) is regular in a strip ¢, < Re ¢ < 0, whereas 

3ergman requires only regularity in a partial neighbourhood of ¢ = 0. Against this must 
be set the fact that Bergman’s form is established only when @ is restricted as in (14), 
while in the Chaplygin form @ is unrestricted. 

*See particularly [5], p. 258, un-numbered equation following (4); here y* is defined in (1.6), where 


z-1/2 ig the same as V(r) of the present paper. 
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THE DESIGN OF TWO-DIMENSIONAL CONTRACTION SECTIONS* 
By PAUL A. LIBBY anp HOWARD R. REISS (Polytechnic Institute of Brooklyn) 


1. Introduction. The design of the expansion section of two-dimensional or axially 
symmetric supersonic nozzles has been the subject of considerable investigation. The 
methods to be used in their design are widely understood and readily applied. These 
are all based on the assumption of uniform, locally sonic flow at the minimum section. 
However, the literature on the design of a contraction section which will produce this 
uniform flow appears to be limited. 

Tsien [1], Szezeniowski [2] and Smith and Wang [3] have presented methods for the 
design of axially symmetric contraction cones for an inviscid, incompressible fluid. There 
are, however, cases in which the contraction section is essentially two-dimensional and in 
which the inlet flow is approximately axial and uniform. For such contraction sections 
the authors were able to find no literature, although the above referenced methods 
could be applied to this case. However, the hodograph method [4] for incompressible, 
two-dimensional potential flow appears to offer a simpler solution to this problem. 

2. Development. In this method the actual contraction section is taken from an 
infinite channel with asymptotes of y = +a at x — —o and of y = +8 at r >. 
The velocity at these extremes is in the x direction and is equal to a and b respectively. 
It is required to find a contour connecting these asymptotes in such a fashion that the 
resultant velocity is everywhere monotonically increasing in the direction of flow. This 
requirement assures for both the incompressible and compressible flow that boundary 
layer separation and local compressibility effects will be avoided (cf. reference 1). 

Although there may be chosen an infinite variety of functions, which give streamlines 
connecting points (a, 0) and (b, 0) in the hodograph (u, v) plane and which satisfy this 
requirement, a simple set may be found as follows: Consider the complex potential 
F=¢o+ w 


F(w) = c[In (W — 6b) — In (W — a)], (1) 


where @ is the usual complex conjugate velocity, uw — iv, and where c is a positive, real 
constant. Considering the real and imaginary parts of Eq. (1), one finds that 


o/c = (1/2) In {[(u — a)? + v7 J/[(u — b)? + v' J}, (2a) 


W/e = tan™’ [—v/(u — a)] — tan™’ [—v/(u — D)]. (2b) 
It is of interest to examine the streamlines in the W plane. From Eq. (2b) one obtains 
the equation 
v/c ad 6, = 6, ait 


where these angles are shown in Fig. 1. An examination of this figure indicates that the 
streamlines are ares of circles passing through the points (a, 0) and (6, 0). Since the 
velocity at any point on such a streamline is represented in the W plane by a vector to 
that point from the origin, only those streamlines corresponding to 7/2 < | ¥/c| < « 
will yield wall shapes in the z plane with monotonically increasing velocities. The arcs 


*Received June 5, 1950. 
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corresponding to | ¥/c | < 2/2 yield, in the z plane, streamlines along which the velocity 
is locally decreasing in two regions. Figure 1 also indicates that the centerline of sym- 


-V 














Fic. 1. Streamlines in the w plane. 


metry of the contraction section corresponds to ¥/c = +x. The x and y coordinates of 
the streamlines will now be determined. 
Solving for w from Eq. (1), one obtains 
m= [b — aexp(—F/o)]/[1 — exp (—F/o)]. (3) 
From the definition of F, the complex variable z = x + iy is given by the transformation 
equation 


c= [| @ "dF +C,, (4) 
or from Eq. (3) 
z= / {{1 — exp (—F/c)]/[b — a exp (—F/o)]} dF +C,, (5) 
which, upon integration gives 
2 = (F/b) — (c/b)[(b/a) — 1] In [(6/a) — exp (—F/)] + Ci, (6) 





where C, and Cj are arbitrary complex constants of integration. 
d tn) 
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Now with F = ¢ + w@ one obtains, by equating the real and imaginary parts of 
each side of Eq. (6), 


x = (¢/b) — (c/b)[(b/a) — 1] In {1 + [a exp (—¢/c)/b}’ 
(7a) 


— [2.a cos (¥/e) exp (—#/e)/b]}'* + Cs , 


y= (y/b) —_ (c 'b)[(b/a) —_ 1] 
(7b) 


- tan” {[exp (—¢/c) sin (¥/c)]/[(b/a) — exp (—¢/c) cos ¥/c]} + C: , 
where again C, and C, are arbitrary real constants. To obtain the streamlines of the 
flow given by Eqs. (7a) and (7b) one may consider them, for ¥ = constant, two parametric 


equations giving x and y in terms of ¢. 
To determine c and C, , Eqs. (7a) and (7b) are examined. For¢ — —”, 1 —@ 


provided b/a > 1, and 
y — (p/b) — (c/b)[(b/a) — 1]l+a — (¥/e)] + C2, (8a) 
where the upper and lower signs correspond to ¥/c 2 0 respectively. Furthermore, for 
ro) — Oo .F —o, 
i (y/b) + C, ’ (8b) 
where the arctangent is chosen in the first and second quadrant for ¥/e > O and in 


the third and fourth quadrant for ¥/c < 0. 
The constant C, is determined so that ¥/c = +7 along the z axis (y = 0). Thus 


from Eq. (7b) 


C, = ¥Fer/b, (9a) 
giving as ¢—>—-2, y—[(¥/c) F r)/(a/c), (9b) 
and as oo, y — [(v/c) = x) /(b/c). (9c) 


Now if it is required that as 7 > —© (6 — —~©&), y — Fa for ¥/c = +p, i.e., 
that the contraction section walls correspond to ¥/c = +y> , then from Eq. (9b) 


a = (—% + 2)/(a/c) 


or 

c = aa/(—% + =). (10a) 
Then as x ~& (¢ >), from Eqs. (9c) and (10a) it follows that 

y — Faa/b = F8£. (10b) 
Thus if = (a/b) exp (—¢/c), Eqs. (7a) and (7b) with Eqs. (9a), (10a) and (10b) yield 
for the wall contour in the upper half of the z plane (¥/c = —y») denoted by the subscript 
zero 
£, = (2)/a) 

(11a) 


= —[(B/a)/(—vo + m]{In & + 1/2[(@/8) — 1] n(1 + & — 26 cos Y)} 
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=) 
Go 


and 
no = (Yo/a) = [(B/a)/(— Yo + 7)] 
(11b) 
{—yY. — [(a/8) — 1] tan’ [—@sin y,/(1 — © cos y%)] + 7}. 


The arbitrary constant C, in Eq. (7a) leads only to a translation of the streamlines along 
the x or £ axis and thus may conveniently be set equal to zero. 

Equations (lla) and (11b) give the final contour in the upper half plane in terms of 
the parameter (0 < 6 <o~). Note that & — —o and m» — 1 as @—~o, and &, ~@ 
and  — 8/a as ® — 0. In order to obtain the streamlines inside the wall contour, ¥/c 
may be permitted to vary from its wall value, —y , to —7. It might be pointed out that 
for asymmetrical channels y) would simply take on different values in the upper and 
lower halves of the z plane. 

3. Numerical Example. In Fig. 2 the wall contours for a/8 = 4 with y = 7/2, 
27/3 and 7/3, are shown. Note that the latter contour, for which the velocity is not 





—-— W=1/3 
W = 1/2 
alt ati a Wo = 2TI/3 














“3 “— =| e) | 
$0 


Fic. 2. Wall contours. Contraction ratio = 4. 


monotonically increasing, appears to be reasonable and might well be drawn and selected 
for a contraction cone if no analysis were carried out. It will be noted that the asymptotic 
values, 1 and 8/a are approached quickly. Criteria of closeness to the asymptotes, of 
uniformity of the u velocity, or of the smallness of the v velocity can be applied to estab- 
lish the values of & determining the length of the section. 
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A NOTE ON BATEMAN’S VARIATIONAL PRINCIPLE 
FOR COMPRESSIBLE FLUID FLOW* 


By CHI-TEH WANG (New York University) 


1. Introduction. The variational principle for a compressible fluid was first studied by 
Hargreaves [1], who showed that the integrand of the variational integral is a linear 
function of the pressure. A variational principle for an inviscid compressible fluid was 
formulated by Bateman [2]. A study of Bateman’s work, however, shows that his varia- 
tional principle is applicable only when the domain of the flow is finite. A large class 
of aerodynamic problems require the study of a flow field which extends to infinity. In 
such cases, Bateman’s principle must be modified. This fact has already been noted in 
references [3] and [4], in which the Rayleigh-Ritz method was used in the approximate 
solution of compressible flows past arbitrary bodies. The formulation of a suitable varia- 
tional principle in these references was however carried out in connection with the 
particular problems considered, so that the derivation appears to be in a rather restricted 
form. In this note, a more general formulation is presented and the resulting variational 
integral is written in a more general form. The author is indebted to Professor K. O. 
Friedrichs for his kind suggestions and discussions. 

2. Bateman’s variational principle. For steady, inviscid, irrotational compressible 
flow, the governing differential equation is 


E (5 | 3 _ 9b 8 8 _ (1) 


Ox; Ox; OX; OX; OX; OX; OX; 7 
W he e 
2 = 1 2 Od 0d ° . 
= <a (<f, ond az. ae.) (i, J) k= 1, 2, 3). (2) 
“ “Hk “k 


Cartesian coordinates, y is the ratio of specific heats, q, is the maximum attainable 
velocity in the flow. A repetition of the subscripts in the above expression indicates 


In Eqs. (1) and (2), a is the velocity of sound, ¢ is the velocity potential, x; are the 
summation. 
Bateman’s problem is to show that the variational integral 
1, = [ p@av (3) 
“V 
has Eq. (1) as its Euler’s equation, where p is the pressure, d V is the elementary volume, 


and the integration is extended to the whole volume of the fluid. 
For barometric fluid, p may be written as 


p= A+ Bp’, 
where p is the density, and A, B are constants. In terms of ¢, one obtains 
> _ db dd ) 
= / 1g? —- —— S 4 
tail o(¢: Ax; ax, (4) 
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where 
C= BMO-M'y — 1/2y)* and k = y/(y — 1). 
With the expression of p as given by (4), the first variation of J, may be carried out 
with respect to ¢ and the condition 6J, = 0 leads to 
"8 ( a ) em ag , 
—~{g——) adv — dS = 0 é 
| ” Ox (, Ox; ov + ] ds on as : (5) 


i “S,+S8s, 


where 6 indicates first variation, 0/dn is the derivative in the inward normal direction, 
S, + S, is the surface that encloses the volume V, S, denotes the stream surface and 
S, denotes the boundary surface at infinity. Since 6¢ is arbitrary in V, the condition 
6I, = 0 gives the continuity equation 


oO os _ 
ax, (, ax,) 0 


as the Euler’s equation of the variational integral (3). 
If the domain is finite and the boundary surfaces are stream surfaces, no condition 


has to be imposed on ¢ and d¢/dn = 0 follows as the natural boundary condition. If 
the domain is infinite, 0¢/d8n = O on S, can still be concluded from the condition 


fs, 66 p(d¢/dn) dS = 0. At infinity @ must be prescribed, and since S, is an infinite 
surface, the vanishing of 67, requires that ¢ must be prescribed to an order of magnitude 
so that [s, 66 p(d¢/dn) ds = 0. This however is not the case in fluid dynamics problems. 
To clarify this point, let us consider the two-dimensional case. At infinity, the ad- 
mitted velocity potential is required to behave as follows 
; r + lana -1 -1 0° —2 
@ = Ur cos 6 — on 6+ (U+ A,)r- cos 6+ A,r sin 0+ O(r”), (6) 


where 7, 6 are the polar coordinates, K is the strength of circulation A, and A, are to 


be determined and O(r~*) represents terms of the order 1/r° or higher. Thus 

6¢ = 6A,r™' cos 6 + 6A,r' sin 6+ O(r”), 
where 6A, and 6A, are arbitrary. Writing p and d¢/dn in terms of ¢ as given in (6) and 
integrating, one obtains 


nG=a2"r 


[otter | 


| i@ pt r ao} = —poUrbA, , (7) 
where p, is the density at infinity, s, is the boundary curve and ds is the elementary 
length. The vanishing of 6/, then requires that p) U must be zero because 6A, is arbitrary. 
This however is not possible. It is therefore clear that Bateman’s variational principle 
is not applicable to flows in which the domain extends to infinity. The appropriate 


principle in this case should be 
i | p(o) dS + rd, | = 0), (8) 


where S is the surface of the domain and dS is the elementary surface. 
3. A variational principle for steady, irrotational compressible flow with infinite do- 
main. In applying the Rayleigh-Ritz method to the approximate solution of compressible 
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flow problems, it is found that for most problems the velocity potential @ may be written 


in the form 
ty) — Pi: + deo ? (9) 


where ¢, denotes the velocity potential of the corresponding incompressible flow and 
¢, denotes the remaining part due to the compressibility effect. Substituting ¢ as given 
by (9) into Eq. (4) and carrying out the expansion, the expression for p may be written 
in the following form 


— Y dg, 0g, fo)" _ on Og. re 
= 44 ee - nad “ee” (10) 
Since ¢, is a definite function, 5¢, is zero, and hence 
(2 _ oy a1)" —_ 
6 [ | 4 + of ax, ax; dv = 0, (11) 
, 0d; 0d» “awed 2th _ : 
. I, " O22; OX; adil i bo. Po 9 dS [ La” dx, 2, dV 


(12) 


dg, 
~f 62 Po ry dS. 
The last step in (12) is obtained because ¢, satisfied the Laplace equation 0°¢,/dz7,0x; = 
0 and on the stream surface S, , d¢,/d8n = 0 and thus fs, 6¢2(0¢,/dn) dS = 0. 

As long as the Euler equation of a variational integral is not affected, it is permissible 
to change the original integral by adding or subtracting other integrals. Since (11) is 
zero and (12) gives only a boundary integral, the following variational integral will have 
the same Euler’s equation as Bateman’s integral (4). 


it j f Oe 7 ig — dg, 261)" ]} |p 0g; Og. 
= } \Pe E t o(«: dx; Ox; eV + | ae, dn,°" a 
Noting that [A + C(qi, — (0¢,;/d2x;)(0¢,/dx,)"] = p(¢,) and integrating the second 
integral in (13) by Green’s formula, (13) becomes 


= | w@ - poniav — | dm Sas. (14) 


In the above expression, p(¢,) subtracted thusly insures the boundedness of J, . The 
last integral in (14) must be subtracted so that the boundary integral vanishes om 
the first variation of (14) is taken. In the two-dimensional case, the first variation of 
the last integral in (14) indeed reduces to (7). 

In references [2]-[9], variational methods have been carried out to solve compressible 
flow problems following the Rayleigh-Ritz, the Galerkin, and Biezieno-Koch procedures. 
In all the problems solved, excellent results were obtained. In the case of the Galerkin 
method and the Biezieno-Koch method, the formulation of a variational principle is not 
necessary. However, in performing the numerical computation for potential flows past 
arbitrary bodies, it was found that the Rayleigh-Ritz method requires the least amount 


of labor. 
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ON A GEOMETRICAL METHOD OF DERIVING THREE-DIMENSIONAL 

HARMONIC FLOWS FROM TWO-DIMENSIONAL ONES* 
By AUREL WINTNER (The Johns Hopkins University) 

The Flow Operators 2. Let D be a domain in a (u,v)-plane, E a domain in an (2,y,2)- 

space, and let g = ¢(u, v), ¥ = ¥(z, y, z) denote (real-valued and regular) solutions of 


Ao.w = 0, Aww = 0 on D, E, respectively, where A, and A; denote the two- and three- 
dimensional euclidean Laplace operators, 


92 /4.,2 42 2 2/4,.2 42 2 2 2 
A, = 0°/du’ + d°/av and A; = 0°/dx + d°/dy + 0°/dz. (1) 

All harmonic functions ¢ are accessible in principle, since all of them are given by 
g(u, v) = Rex(w), where x is any function which is regular-analytic in w = u + w on 
D. In contrast, there does not exist anything like this rule for the harmonic functions 
¥(x, y, 2) on a three-dimensional Z. Hence it is natural to ask for flow operators, say 


2 = 2(D), which, from every regular solution g = ¢g(u, v) of A,g = 0 on a two-dimen- 
sional (u,v)-domain D, will manufacture a regular solution, 

V(x, Y; 2) = Qu, Vv), (2) 
of A,wwW = 0 on a three-dimensional (x,y,z)-domain E = E(D). The latter should not 


depend on the particular choice of the function g(u, v), but merely on the operator 
Q = 2(D) and on the domain D on which ¢(u, v) is supposed to be harmonic. 

A trivial instance of such “harmonic flow operators” @ is supplied by the cylindrical 
flow which, from a given ¢(u, v), manufactures the corresponding y(z, y, 2) as follows: 


V(x, y, 2) = of, y). (3) 
In fact, (3) is of the type (2), since A;y(z, y, z) = A.w(za, y) if ay /dz’ = 0; cf. (1). The 


*Received June 30, 1950. 
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E = E(D) belonging to the case (3) of (2) is the interior of the infinite cylinder having 
the (orthogonal) cross-section D. 

The History of the Problem. It turns out that this trivial flow operator Q, that defined 
by (3), is only the simplest instance of a large class of flow operators 2 which have the 
property specified in connection with (2). In fact, the general class of harmonic flow 
operators 2 was discovered by Weingarten sixty years ago.’ His results are stated, 
however, in the language of the differential geometry of minimum surfaces, and seem to 
have remained unnoticed (except, possibly,” by F. Klein). In any case, the principal 
result of Weingarten was rediscovered by Levi-Civita a decade later,* by using, instead 
of the theory of minimum surfaces, extremely heavy explicit calculations, which he 
only was able to carry out due to Ricci’s and his absolute differential calculus. In addi- 
tion, these results of Levi-Civita are in a section which is not along the main lines of 
his long paper. Finally, the latter was published in a periodical not available in most 
libraries. 

These circumstances notwithstanding, it is somewhat surprising that the Q-flows in 
question seem to have been forgotten for half a century. On the other hand, the nature 
and the generality of the result appear to be elastic enough to admit its explicit use in 
solving certain three-dimensional boundary-value problems which could hardly be 
attacked in any other way. Under these circumstances, and because only the proof, 
but not the final result, is quite involved, the criterion in question seems to be worth 
recording here in a form which, it is hoped, can be followed by those who apply mathe- 
matics. 

(i) Line Congruences. Let T be a 2-parameter family of straight lines which, with 
reference to a three-dimensional (x,y,z,)-domain EF, have the following property: If P 
is any point of Z, then T contains exactly one straight line, say y = yp , through P. In 
the terminology of the differential geometry of “rays” y, such a collection T is called a 
“line congruence” (on £). It will be assumed that that portion of each of the straight 
lines which is contained in £ is a connected set (i.e., that it is a segment or a half-line, 
possibly the entire line y). 

Analytica!!y, every sufficiently small portion of I can be described by first choosing 
inside E a surface, say S, and then introducing on S Gaussian parameters, say u and », 
in an arbitrary way, subject only to the restriction that there be a one-to-one (and 
sufficiently differentiable) correspondence between the points of S and the points of a 
domain, say D, in a euclidean (u,v)-plane. In fact, if S and its (u,v)-parametrization 
are fixed, then a straight line y belonging to T can be individualized by placing y = 
y(u, v), where u, v are the parameter values of that point of S from which y is issued. 

(ii) Rectilinear Flows. In this manner, the line congruence I defines a rectilinear 
flow from the surface S into the three-dimensional (z,y,z)-domain FE (which contains S). 


1J. Weingarten, Ueber particulire Integrale der Differentialgleichung AV = 0 und eine mit der Theo- 
rie der Minimalfliichen zusammenhiingende Gattung von Fliissigkeitsbewegungen, Nachrichten von der 
Kéniglichen Gesellschaft der Wissenschaften zu Géttingen, 1890, pp. 313-335. 

2Cf. the footnote on p. 117 of Levi-Civita’s paper, referred to in the next footnote. The common 
element in the treatments of Weingarten and Levi-Civita is Jacobi’s theorem concerning the (complex- 
valued) solutions f = f(z,y,z) of (af/daxr)? + (af/dy)? + (af/daz)? = 0, a theorem to which Levi-Civita’s 
attention was called by Klein (ef. the footnote mentioned before). 

8T. Levi-Civita, Tipi di potenziali che si possono far dipendere da due sole coordinate, Memorie della 
teale Accademia delle Scienze de Torino, ser. 2, vol. 49 (1900), pp. 105-152; more particularly, pp. 138-139. 
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Correspondingly, if ¢ is any function of the position (wu, v) on a two-dimensional domain, 
D, in the euclidean (u,v)-plane, with u, v as Gaussian parameters on S, then a function, 
say y, of the position (x, y, z) on E will be defined by placing 


V(x, y, 2) = ou, v) if (2, y, z) is on y(u, v). (4) 


It is understood that y(u, v) in (4) denotes that straight line, y, of ! which reaches 
the surface S at the point having the Gaussian coordinates (u, v). 

(iii) Laplacean Rectilinear Flows. Clearly, the operation (4) is of the type (2). But 
it was not assumed in (4) that the given function, ¢(u, v), satisfies A.g = 0. Moreover, 
if A.g = 0 is assumed, it is not in general true that the function y(z, y, z) defined by (4) 
will satisfy A,y = 0. The operators A, , A; are those defined by (1). 

This problem is solved by the result of Weingarten and Levi-Civita, referred to 
above. It states that, if g(u, v) is any (continuous) solution of A,g = 0 on a (u,v)-domain 
D, then the function ¥(2, y, z) defined by (4) will be a solution of A;y = 0 on the corre- 
sponding (z,y,z)-domain E = E>, provided that the rectilinear flow, on which the assign- 
ment (4) is based, is derived from a line congruence I which is an isotropic line con- 
gruence. 

The latter notion, to be defined under (iv) below, depends only on the collection of 
lines which constitute I’, rather than on the particular choice of the surface S and of the 
Gaussian parametrization (u, v) of S. Consequently, neither of these choices matters, 
even though both of them occur in (4). 

(iv) Isotropy. There remains to be defined the notion of an isotropic line congruence, 
mentioned above. The customary definition of this notion, introduced by Ribacour,* 
is quite involved. In what follows, it will be replaced by another definition, one which 
can be stated more easily and which, as a matter of proof,” turns out to be equivalent 
to the usual definition. 

Consider in the (z,y,z)-space two (sufficiently small and sufficiently differentiable) 
pieces of surfaces, say S and T, and let P — Q = Q» be a one-to-one (and sufficiently 
differentiable) mapping of the points, P, of S on the points, Q, of T. Suppose that every 
S-point, P, is distinct from its T-image, Q = Qp , i.e., that every pair P, Q determines 
a straight line y = y (P; Q) and therefore the entire mapping, S — T, defines a certain 
line congruence, say I’. Suppose further that the mapping S — T' defining this T has 
the following two properties: 

(a) P—Q = Q> is an isometric mapping of S on T' (as to this notion, cf. (v) below) 
and 

(8) | PQ| = const., i.e., the euclidean length of the segment joining P with Q = Qp 
(in the (x,y,z)-space) is independent of the choice of P on S. 

If these conditions are satisfied by the mapping S — T’, then the line congruence I, 
mentioned before (a), is an isotropic line congruence. Conversely, every isotropic line 
congruence can be obtained in this manner, that is, by suitably choosing two surfaces, 
S and T, and a mapping, S — 7’, which satisfies both of the above conditions, (a) and (8). 

(v) Isometry. For the sake of completeness, the notion of isometry, occurring above 
in (a), remains to be explained. This notion is familiar from the elements of the theory 
of surfaces. It is defined as follows: 


4L. Bianchi, Lezioni di geometria differenziale, 2nd edition, vol. 1 (1902), p. 302. 
5Loc. cit.,4 p. 304 
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A one-to-one mapping of the points, P, of the surface, S, on the points, Q = Qp, 
of another surface, 7’, is called isometric if, by virtue of the mapping, the squared line- 
element, ds”, on S (i.e., the so-called “metric” or “first fundamental form”, g;, du’ du‘, 
on S) coincides with the squared line-element on T. 

In order that this be the case, it is necessary that the Gaussian curvature of S at 
any point P be identical with the Gaussian curvature of 7 at the corresponding point 
Q = Q> . This necessary condition is sufficient as well in the particular case of surfaces 
of constant Gaussian curvature. Since cylinders are of constant Gaussian curvature, the 
constant being 0, this embeds the trivial example (3) into the general theory. 


FREE LONGITUDINAL VIBRATION OF A PROLATE ELLIPSOID, 
CLAMPED CENTRALLY* 


By JAMES S. KOUVELITES (Sloane Physics Laboratory, Yale University**) 


Introduction. It was originally shown by Poisson’ that the magnetization is uniform 
and parallel throughout the interior of an ellipsoid of magnetically isotropic material 
placed in a previously uniform and parallel static field. Since this unique property of 
the ellipsoid allows an exact calculation of its demagnetizing factor as well as a simple 
field analysis, the prolate spheroid was chosen as the appropriate shape for the speci- 
mens whose magnetostrictive vibrational properties are being studied in this laboratory.” 
Although the magnetostrictive vibration is forced rather than free, the damping factor 
has been found to be so small for the materials tested, that through an investigation 
of the free vibration, a reasonably good approximation to some of the resonance phe- 
nomena may be obtained. 

Analysis and discussion. For later comparison with the corresponding expressions 
for the ellipsoid, the differential equation of motion, its integrated solution, and a 
subsidiary expression for the frequencies of resonance for free longitudinal vibration of 
a bar of constant cross-section, clamped centrally,’ are 


aE _ pdt 
Ox’ —_— k or (1) 
1/2 
& = Asin (.(2) r) sin (wt + ¢) (2) 
, n k 1/2 " 
Sr, = Z\7 ’ n = 1, 3,5, -*> (3) 
4a \p 


“Received Sept. 26, 1950. 

** Assisted by the Office of Naval Research. 

tAs described in J. C. Maxwell, A treatise on electricity and magnetism, vol. II, 3rd ed., Clarendon 
Press, Oxford, 1892, pp. 66-69. 

2J. S. Kouvelites and L. W. McKeehan, Magnetostrictive vibration of prolate spheroids. Preliminary 
measurements, Rev. Sci. Instr., in press. 

’Rayleigh, Theory of sound, vol. I, §150. 
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where x is the coordinate along the longitudinal axis (the origin being at the center of 
the bar), ¢ the time, & the longitudinal displacement, p the volume density, k the co- 
efficient of elasticity, a the half-length of the bar, w the angular frequency, A and ¢ 
arbitrary amplitude and phase constants. 

For an elastically-isotropic prolate ellipsoid with two eccentricities nearly equal to 
unity as in the physical problem of magnetostrictive vibration, it appears reasonable 
to assume that the equiphase surfaces of £ are planes perpendicular to the major axis.* 
Applying the second law of motion to a transverse plane lamina of thickness dz of an 


ellipsoid bounded by the surface 


zx y Zz — : - 
+ 3+3=1, where 2a is the major axis, 


oF at x 
— dx = —3 dm, where dm = ryzp dx = rie 1 — ay dx. 
Ox dt E a 


Since Hooke’s law may be stated as 
a°\, d€ 
F = re 1 = a — ; 


the equation of motion becomes 


2 
6 


Oo} 2x 0& poet (4) 
Ox" a—x ox kot’ 


Str 


which is a definite departure from the familiar wave equation (1) appropriate to a 
channel of uniform cross section. 

The separation of the variables, x and t, is accomplished as usual by letting § = 
X(zx)T(t), which is substituted into equation (4). Both sides of the resultant expression 


may then be set equal to the constant, —w p/k, yielding: 
d°T 5 
—— twill =0 (5) 
dt 
dX 2x dx ~~ 
—s 2 + wo ‘ X = 0 (6) 
dx a—x dx k 
Thus, & = X sin (wt + ¢) where X may be obtained as an explicit function of x by in- 


serting it into Eq. (6) in the form of the power series, 
p p+ \pts 
- HY (x oi —_ 
xX A?) + 4, } go> + 442) deans (7) 
\a \d, a 
As a result, it is found that, if p = 0, both A, and A, have arbitrary, finite values and 
the coefficients of Eq. (7) satisfy the condition 


8 N 
— ——_ _ — _ ~ ——— ( 9 ——_ Zz 
Aus = GA. +7  (A,-2 — Ai), (8) 


where N = w’a’p/k 





‘See Rayleigh’, Vol. II, §265 and P. M. Morse, Vibration and sound, McGraw-Hill Book Co., New 
York, 1936, pp. 216-217 for analogous treatment of sound transmission in tubes of varying cross-section. 
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Moreover, since the physical situation requires that the ellipsoid be clamped centrally 
(i.e. £ = O at x = O at any time), Ay , A. , A, , etc., must be zero. The series solution 
for X, subject_to Eq. (8), then becomes 


X= A,(£) + a,(2) + a2) fee, (9) 


0 0.2 0.4 0.6 0.8 1.0 


+0.8 T ] T T SF i -, y +0.8 
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Fia. 1. 


There is no force exerted on the ends of the ellipsoid. Consequently, the second 
boundary condition of the problem is that dX/dz = 0 at x = +a. The particular values 
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of N which satisfy this condition, are found by a trial and error process from Eqs. (8) 
and (9). It is interesting to note that Eq. (9), for x = ta, is divergent for all values 
of N except those corresponding to the modes of vibration satisfying the second boundary 
condition. 

In Fig. 1, X/A, , has been plotted as a function of x/a for various values of N, the 
thickly-drawn curves representing the “standing-wave” patterns for the fundamental 
mode and its first two overtones. It was sufficient to show the deflection variations over 
only one half the length of the ellipsoid because there is symmetry about the plane 
x= 0. 

The frequencies of resonance for the ellipsoid may thus be expressed as 


1/2 
fe == (*) where n & 1.38, 3.49, 5.33, ---, (10) 
4a \p 

indicating the inharmonic nature of the overtones. However, a comparison between 
Eqs. (3) and (10) reveals that the consecutive overtones of the vibrating ellipsoid have 
ratios closer and closer to the ratios of adjacent odd integers as the value of n increases. 

Acknowledgements. I wish to thank Professor L. W. McKeehan of this laboratory, 
whom I am assisting in the study of magnetostrictive vibration, and Professor F. J. 
Beck, Jr., of the Dunham Laboratory of Electrical Engineering, for encouraging discus- 


sions of this problem. 


THE LEAST SQUARES SOLUTION FOR A SET OF 
COMPLEX LINEAR EQUATIONS* 


By R. TURETSKY (Aberdeen Proving Ground) 
Consider the set of m observational equations whose matrix representation is 


Az~ wv, (1) 


where A is an m X n (m > n) matrix of rank n whose elements are prescribed complex 
quantities, while w is an m X 1 matrix. We seek that z (an n X 1 matrix) which mini- 
mizes the sum of the squares of the absolute values of the components of the vector 
Az — w. 

SettA = B+ 70,2 =2+ ty, w = ut w, where B, C, 2, y, u, and v are real matrices. 


Then Eq. (1) is equivalent to 


[> 


~ 2) 


Cc B/\y v 


This is the matrix representation for a set of real observational equations. To obtain 
the normal equations,’ we multiply on the left by the transpose of the coefficient matrix 


*Received July 28, 1950. 
1See Whittaker and Robinson, Calculus of observations, Chapter IX. 
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of the unknowns. Letting primes denote the transpose, there results 


P—Q\/z B’ CO Vse 
= , (3) 
Q P/\y —C’ B’'/ Nov 
where 
P = B’B + C'C 
and (4) 


Q = C’B — BC 


It will be noted that P is symmetric while Q is skew-symmetric, which greatly facilitates 
the computation. 
The solution of Eq. (3) is 


x —Q\"'/B’ C’'\fu 
y Q P C’ BI Nu 
Ordinarily the inverse of the 2n X 2n matrix 
P -Q 
Q P 
recourse may be had to its special structure to note that 


is obtained directly. However, 


the inverse is of the form 


where 


R = (P + QP"'Q)" 
(6) 
S = P"QR 


The computation is thus simplified, since we have to compute the inverses of two n X n 
matrices rather than that of one 2n X 2n matrix. Moreover, advantage can be taken 
of the fact that R is symmetric while S is skew-symmetric. 

Finally, the multiplication of 


R S 
—-S R 

by 
B Cl 


—-C’ B’ 
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required on the right hand side of Eq. (5) leads to the matrix 
T -U 
U g 
where 
T = RB’ — sC’ 
UV 


Multiplication by a 2m X 2n matrix is thus required only in the finz! step. 


Il 


RC’ + sB’ 





CLARITY is the outstanding advantage of... 
MATHEMATICS of RELATIVITY 


By G. Y. Ratnicu, University of Michigan. Professor Rainich takes the well-known facts of relativity and 
arranges them so that the theory appears as a consistent whole. To make his presentation as clear as possible, 
the author uses the following approach: 

He stresses the mathematical aspects of relativity, beginning with a discussion of the inverse square law 
in terms of simple calculus. The more complicated situations are introduced step by step to bring in the more 
sophisticated mathematical tools. 

He takes up the difficult points (such as changes in the fundamental concepts and refinements of the mathe- 
matical techniques) only after the need for them is made clear. In this way he makes the book easier to 
follow without sacrificing rigor. 


Mathematics of Relativity goes beyond earlier books in one important respect: In section 19 Professor 


Rainich demonstrates that both gravitation and electromagnetism fit exactly into the original theory of curved 
{ volume in the Applied Mathematics Series, edited by 1. S. Sokolnikoff. 1950. 173 pages. $3.50. 
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space. 


By Henry L. Lancuaar, University of Illinois. A systematic, 


DIMENSIONAL ANALYSIS 
thorough treatment of the principles of dimensional analysis, a 
and THEORY of MODELS subject the author defines as “a method by which we deduce 


information about a phenomenon from the single premise that 
the phenomenon can be described as a dimensionally correct 
equation among certain variables.” 

Dr. Langhaar stresses the engineering applications of dimen- 
sional analysis and the engineering use of models. The material 
is organized so that the engineer can apply the theory without 
understanding the mathematics completely. April 1951. 

100 pages. $4.00. 


By L. H. C. Tippett, British Cotton Industry Research Assocta- 


TECHNOLOGICAL APPLICATIONS 
tion. Stresses the practical points that come up in applying 
of STATISTICS statistics to technology, and the practical importance of the 


mathematical assumptions involved. Few men are as well quali- 
fied as Mr. Tippett to write on this subject. He is an inter- 
nationally-known authority on the applications of statistics in 
industry, and particularly in industrial research. 

By presenting his material directly, without repetition, Mr. 
Tippett succeeds in covering control charts, sampling inspection, 
analysis of variance, and design of experimentation—all in a 
comparatively short book. One of the Wiley Publications in 


Statistics, Walter A. Shewhart, Editor. 1950. 189 pages. $3.50. 
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